#COMPILE DLL "sdsun.dll"

' This file is part of the SunDLL.  The SunDLL is free software: you can
' redistribute it and/or modify it under the terms of the GNU General Public
' License as published by the Free Software Foundation, version 2.
'
' This program is distributed in the hope that it will be useful, but WITHOUT
' ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
' FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
' details.
'
' You should have received a copy of the GNU General Public License along with
' this program; if not, write to the Free Software Foundation, Inc., 51
' Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
'
' Copyright David R. Pratten


' 03.01.02
'   Updated DT-UT FOR 2002-2006 FROM http://hpiers.obspm.fr/eoppc/bul/bulc/UTC-TAI.history
' 03.01.01
'   First Freeware version

' WARNING WARNING WARNING
'
'do NOT compare for  x = %TRUE directly as %TRUE = 1 and PB returns -1 for true result
'  USE 'ISTRUE x'
'do not say 'if not x'
'   say 'if isfalse x' instead
' WARNING WARNING WARNING
'

'
' Comment this line to make DLL otherwise it is an exe that could be debugged
'%DEBUG = 1
'


#REGISTER NONE
#IF %DEF(%DEBUG)
    #DEBUG ERROR ON
#ENDIF


OPTION EXPLICIT

#INCLUDE "WIN32API.INC"
#INCLUDE "registry.inc"
#INCLUDE "md5.inc"
#INCLUDE "crypto.inc"
#INCLUDE "commequa.inc"
'
' Equates
'

GLOBAL EarliestMJD AS DOUBLE
GLOBAL LatestMJD AS DOUBLE


GLOBAL SDYEOTMisJoin AS EXT
GLOBAL SDYDeclMisJoin AS EXT

' Globals for security
GLOBAL Serial AS STRING
GLOBAL KeyFileLevel AS LONG
GLOBAL KeyType AS STRING
GLOBAL Company AS STRING
GLOBAL Email AS STRING
GLOBAL TimeStamp AS STRING
GLOBAL User AS STRING
GLOBAL V() AS DOUBLE




'
' END "equates.inc"
'


DECLARE SUB InitGlobals
DECLARE FUNCTION LoadKeyFile() AS LONG
DECLARE FUNCTION sdInfo() AS LONG



FUNCTION LibMain(BYVAL hInstance AS LONG, _
                 BYVAL Reason    AS LONG, _
                 BYVAL Reserved  AS LONG) EXPORT AS LONG

    '
    ' ENABLE FULL 80bit precision WORK AROUND
    ' Must be first code in module
    '
    #REGISTER NONE
    DIM xyz%
    xyz% = 0
    ASM fstcw xyz%
    ASM OR    xyz%, &B0000001100000000%
    ASM fldcw xyz%


  SELECT CASE Reason

    CASE %DLL_PROCESS_ATTACH


        DIM ret AS LONG
        ret = LoadKeyFile

        InitGlobals

       LibMain = ret   'success?
      EXIT FUNCTION
    CASE %DLL_PROCESS_DETACH
      EXIT FUNCTION
    CASE %DLL_THREAD_ATTACH
      EXIT FUNCTION
    CASE %DLL_THREAD_DETACH
      EXIT FUNCTION

END SELECT
END FUNCTION


#INCLUDE "xml.inc"

#INCLUDE "dtrig.inc"        ' trig library in degrees


#INCLUDE "convert.inc"      ' Unit conversion routines

'#INCLUDE "julian.inc"       ' Calendar routines
FUNCTION CalcJD(BYVAL Year AS LONG, BYVAL  month AS LONG, BYVAL Day AS LONG ) AS EXT
DIM  y AS EXT
DIM m AS EXT
DIM d AS EXT
DIM  A AS EXT
DIM  B AS EXT
DIM  JDpartial AS EXT
DIM xJd AS EXT

IF month < 3 THEN
    m = month + 12
    y = Year - 1
ELSE
    m = month
    y = Year
END IF

A = INT(y / 100)
B = 2 - A + INT(A / 4)
xJD = INT(365.25 * (y + 4716)) + INT(30.6001 * (m + 1)) + Day + B - 1524.5
IF xJD < 2299160.5 THEN
    xJD = xJD - B
END IF
FUNCTION = xJD
END FUNCTION

SUB PartsOfJD(BYVAL JD AS EXT,BYREF YearOfJD AS LONG, BYREF MonthOfJD AS LONG,BYREF DayOfJD AS LONG, BYREF HourOfJD AS LONG, BYREF MinuteOfJD AS LONG, BYREF SecOfJD AS EXT)
DIM Z AS LONG
DIM F AS EXT
DIM jdplus AS EXT

jdplus = jd+.5
Z = INT(jdplus)
F = jdplus - Z

DIM A AS LONG
DIM xalpha AS LONG


IF Z < 2299161 THEN
    A = Z
ELSE
    xalpha = INT((Z-1867216.25)/36524.25)
    A = Z + 1 + xalpha - INT(xalpha/4)
END IF

DIM B AS LONG
DIM C AS LONG
DIM D AS LONG
DIM E AS LONG

B = A  + 1524
C = INT((B-122.1)/365.25)
D = INT(365.25*C)
E = INT((B-D)/30.6001)

DayOfJD = B - D - INT(30.6001*E) 'without part of day

IF E < 14 THEN
    MonthOfJD = E - 1
ELSE
    MonthOfJD = E - 13
END IF

IF MonthOfJD > 2 THEN
   YearOfJD = C - 4716
ELSE
    YearOfJD = C - 4715
END IF

F = F*24
HourOfJD = INT(F)

F = F - HourOfJD

IF F > 0.999999999 THEN
   F = 0
   INCR HourOfJD
END IF

F = F*60
MinuteOfJD = INT(F)

F = F - MinuteOfJD
IF F > .999999999 THEN
   F = 0
   INCR MinuteOfJD
END IF

SecOfJD = F*60 ' only accurate to 5 decimals.

END SUB

FUNCTION JDEMinusJD(BYVAL JD AS EXT) AS EXT
'
' Defines for TD-UT values
'

' UPDATES AT:
' tdut.inc

%tdutstart = 1620
%tdutend = 2006

DIM tdut(%tdutstart:%tdutend) AS STATIC DOUBLE

IF tdut(%tdutstart) = 0 THEN
    #INCLUDE "tdut.inc"
END IF

DIM y AS LONG
DIM m AS LONG
DIM d AS LONG
DIM h AS LONG
DIM mi AS LONG
DIM s AS EXT
PartsOfJD JD,y,mi,d,h,m,s

'MsgBox Str$(JD)+Str$(y)

DIM deltat AS EXT

IF %tdutstart < = y AND y <= %tdutend THEN
   deltat =  tdut(y)/%secinday
   FUNCTION = deltat
   EXIT FUNCTION
END IF
DIM t AS EXT
t = (y-2000)/100


IF y < 948 THEN
   deltat = (2177 + 497*t + 44.1*t^2)/%secinday
'   MsgBox "Deltat"+Str$(deltat*%secinday)
   FUNCTION = deltat
   EXIT FUNCTION
END IF

DIM ErrFactor AS EXT



IF (948 <= y AND y <= 1600) OR %tdutend < y THEN
   deltat = (102 + 102*t + 25.3*t^2)/%secinday
   IF %tdutend < y AND y <= %tdutend+100 THEN
        ' Err Factor is calculated to adjust for matching this formula
        ' To last Observation recorded in "dtut.inc"
        ErrFactor = 40 'seconds
        deltat = deltat + (ErrFactor/100 * (y - %tdutend-100))/%secinday
   END IF

'   MsgBox "Deltat"+Str$(deltat*%secinday)
   FUNCTION = deltat
   EXIT FUNCTION
END IF

IF 1600 < y AND y < %tdutstart THEN
'   MsgBox "Deltat"+Str$(deltat*%secinday)
   deltat = (tdut(%tdutstart) + ( y - %tdutstart))/%secinday
   FUNCTION = deltat
   EXIT FUNCTION
END IF

FUNCTION = %BadParameter ' internal error should never happen

END FUNCTION

FUNCTION JD2JDE(BYVAL JD AS EXT ) AS EXT

FUNCTION = JD + JDEMinusJD(JD)

END FUNCTION

FUNCTION JD2MJD(BYVAL JD AS EXT) AS EXT
FUNCTION = JD - 2400000.5
END FUNCTION

FUNCTION MJD2JD(BYVAL MJD AS EXT) AS EXT
FUNCTION = MJD + 2400000.5

END FUNCTION

FUNCTION LeapYear(BYVAL Year AS LONG) AS LONG
IF Year MOD 4 = 0 THEN
    IF Year > 1582  THEN
            IF Year MOD 100 = 0 AND Year MOD 400 <> 0 THEN
               FUNCTION = %FALSE
               EXIT FUNCTION
            END IF
    END IF
    FUNCTION = %TRUE
    EXIT FUNCTION
END IF
FUNCTION = %FALSE
END FUNCTION


FUNCTION Taum(BYVAL JDe AS EXT ) AS EXT
FUNCTION = (JDe - 2451545) / 365250
END FUNCTION

FUNCTION YMD2SDY (BYVAL Month AS LONG, BYVAL Day AS DOUBLE, BYVAL Hour AS DOUBLE, BYVAL Minute AS DOUBLE, BYVAL Sec AS DOUBLE, BYVAL TimeZone AS DOUBLE) AS EXT
DIM X AS EXT
X = JD2MJD(calcJD(%SDYIndexYear,Month,1))-1.+Day+HMS2Unit(Hour+TimeZone,(Minute),(Sec),%UnitDay)
X = X - JD2MJD(calcJD(%SDYIndexYear-1,12,31))
FUNCTION = Reduce1365(X)
END FUNCTION

FUNCTION YMD2MJD (BYVAL Year AS LONG, BYVAL Month AS LONG, BYVAL  Day AS DOUBLE, BYVAL Hour AS DOUBLE, BYVAL Minute AS DOUBLE, BYVAL Sec AS DOUBLE, BYVAL TimeZone AS DOUBLE) AS DOUBLE 'Ext
DIM X AS EXT
FUNCTION = JD2MJD(calcJD(Year,Month,1))-1.+Day+HMS2Unit(Hour+TimeZone,(Minute),(Sec),%UnitDay)
END FUNCTION

FUNCTION HMS2CTD(BYVAL Hour AS EXT, BYVAL Minute AS EXT, BYVAL Sec AS EXT) AS EXT
DIM x AS EXT
X = HMS2Unit(Hour,Minute,Sec,%UnitDay)
X = X - INT(X)
FUNCTION = X
END FUNCTION

FUNCTION HMST2UTD(BYVAL Hour AS EXT, BYVAL Minute AS EXT, BYVAL Sec AS EXT, TimeZone AS EXT) AS EXT
DIM x AS EXT
X = HMS2Unit(Hour+TimeZone,Minute,Sec,%UnitDay)
X = X - INT(X)
FUNCTION = X
END FUNCTION


FUNCTION SDY2MJD(BYVAL Year AS LONG, BYVAL SDY AS EXT) AS EXT
DIM x AS EXT
x = JD2MJD(calcJD(Year-1,12,31) + SDY)
IF (SDY > 31.0+28.0) AND LeapYear(Year) THEN INCR x

FUNCTION = x
END FUNCTION


'END #INCLUDE "julian.inc"       ' Calendar routines


'#INCLUDE "parmval.inc"      ' parameter validators
FUNCTION BadArcUnit(BYVAL Unit AS LONG) AS LONG
SELECT CASE Unit
    CASE %UnitDegree
        FUNCTION = %FALSE
    CASE %UnitArcMin
        FUNCTION = %FALSE
    CASE %UnitArcSec
        FUNCTION = %FALSE
    CASE %UnitRadian
        FUNCTION = %FALSE
    CASE ELSE
        FUNCTION = %TRUE
END SELECT
END FUNCTION


FUNCTION BadSearchFor(BYVAL SearchFor AS LONG) AS LONG

SELECT CASE SearchFor
    CASE %SearchZeroY
            FUNCTION = %FALSE
    CASE %SearchMaxY
            FUNCTION = %FALSE
    CASE %SearchMinY
            FUNCTION = %FALSE
    CASE ELSE
            FUNCTION = %TRUE
END SELECT
END FUNCTION

FUNCTION BadScale(BYVAL Scale AS LONG) AS LONG
SELECT CASE Scale
    CASE %ScaleMJD
            FUNCTION = %FALSE
    CASE %ScaleSDY
            FUNCTION = %FALSE
    CASE ELSE
            FUNCTION = %TRUE
END SELECT
END FUNCTION

FUNCTION BadFeature(BYVAL Feature AS LONG) AS LONG

SELECT CASE Feature
    CASE %FeatEOTMax1
            FUNCTION = %FALSE
    CASE %FeatEOTMax2
            FUNCTION = %FALSE
 CASE %FeatEOTMin1
            FUNCTION = %FALSE
 CASE %FeatEOTMin2
            FUNCTION = %FALSE
 CASE %FeatEOTZero1
            FUNCTION = %FALSE
 CASE %FeatEOTZero2
            FUNCTION = %FALSE
 CASE %FeatEOTZero3
            FUNCTION = %FALSE
 CASE %FeatEOTZero4
            FUNCTION = %FALSE
 CASE %FeatDeclMax
            FUNCTION = %FALSE
 CASE %FeatDeclMin
            FUNCTION = %FALSE
 CASE %FeatDeclZero1
            FUNCTION = %FALSE
 CASE %FeatDeclZero2
            FUNCTION = %FALSE
 CASE %FeatanalemmaCrossing1
            FUNCTION = %FALSE
 CASE %FeatanalemmaCrossing2
            FUNCTION = %FALSE
 CASE %FeatSolstice1
            FUNCTION = %FALSE
 CASE %FeatSolstice2
            FUNCTION = %FALSE
 CASE %FeatEquinox1
            FUNCTION = %FALSE
 CASE %FeatEquinox2
            FUNCTION = %FALSE
 CASE ELSE
            FUNCTION = %TRUE
END SELECT
END FUNCTION

FUNCTION BadSunType(BYVAL SunType AS LONG) AS LONG
SELECT CASE SunType
    CASE %SunMeanLocal
        FUNCTION = %FALSE
    CASE %SunMeanAtLong
        FUNCTION = %FALSE
    CASE %SunSDYAppAtLong
        FUNCTION = %FALSE
    CASE %SunMJDAppAtLong
        FUNCTION = %FALSE
    CASE ELSE
         FUNCTION = %TRUE
END SELECT
END FUNCTION

FUNCTION BadLongitude(BYVAL Longitude AS DOUBLE) AS LONG
FUNCTION = Longitude <= -180  OR Longitude > 180
END FUNCTION

FUNCTION BadLatitude(BYVAL Latitude AS DOUBLE)  AS LONG
FUNCTION = Latitude < -90. OR 90. < Latitude
END FUNCTION

FUNCTION BadDecl(BYVAL Decl AS DOUBLE)  AS LONG
FUNCTION = Decl < -25. OR 25. < Decl
END FUNCTION

FUNCTION BadMJD(BYVAL MJD AS DOUBLE)  AS LONG
FUNCTION = MJD < EarliestMJD  OR  LatestMJD + 1 <= MJD
END FUNCTION

FUNCTION BadSDY(BYVAL SDY AS DOUBLE)  AS LONG
FUNCTION = SDY < 1  OR  365 + 1 <= SDY
END FUNCTION

FUNCTION BadCTD(BYVAL CTD AS DOUBLE)  AS LONG
FUNCTION = CTD < 0  OR  1 <= CTD
END FUNCTION

FUNCTION BadUTD(BYVAL UTD AS DOUBLE)  AS LONG
FUNCTION = UTD < 0  OR  1 <= UTD
END FUNCTION

FUNCTION BadDate(BYVAL Year AS LONG, BYVAL month AS LONG, BYVAL Day AS DOUBLE)  AS LONG
IF Year < %EarliestYear THEN
   FUNCTION = %TRUE
   EXIT FUNCTION
END IF
IF Year > %LatestYear THEN
   FUNCTION = %TRUE
   EXIT FUNCTION
END IF
IF month < 1 THEN
   FUNCTION = %TRUE
   EXIT FUNCTION
END IF
IF month > 12 THEN
   FUNCTION = %TRUE
   EXIT FUNCTION
END IF
IF Day < 1. THEN
   FUNCTION = %TRUE
   EXIT FUNCTION
END IF
SELECT CASE month
    CASE 1
        IF Day >= 32. THEN
           FUNCTION = %TRUE
           EXIT FUNCTION
        END IF
    CASE 2
        IF LeapYear(Year) THEN
            IF Day >= 30. THEN
               FUNCTION = %TRUE
               EXIT FUNCTION
            END IF
        ELSE
            IF Day >= 29. THEN
               FUNCTION = %TRUE
               EXIT FUNCTION
            END IF
        END IF
    CASE 3
        IF Day >= 32. THEN
           FUNCTION = %TRUE
           EXIT FUNCTION
        END IF
    CASE 4
        IF Day >= 31. THEN
           FUNCTION = %TRUE
           EXIT FUNCTION
        END IF
    CASE 5
        IF Day >= 32. THEN
           FUNCTION = %TRUE
           EXIT FUNCTION
        END IF
    CASE 6
        IF Day >= 31. THEN
           FUNCTION = %TRUE
           EXIT FUNCTION
        END IF
    CASE 7
        IF Day >= 32. THEN
           FUNCTION = %TRUE
           EXIT FUNCTION
        END IF
    CASE 8
        IF Day >= 32. THEN
           FUNCTION = %TRUE
           EXIT FUNCTION
        END IF
    CASE 9
        IF Day >= 31. THEN
           FUNCTION = %TRUE
           EXIT FUNCTION
        END IF
    CASE 10
        IF Day >= 32. THEN
           FUNCTION = %TRUE
           EXIT FUNCTION
        END IF
    CASE 11
        IF Day >= 31. THEN
           FUNCTION = %TRUE
           EXIT FUNCTION
        END IF
    CASE 12
        IF Day >= 32. THEN
           FUNCTION = %TRUE
           EXIT FUNCTION
        END IF

END SELECT
FUNCTION = %FALSE
END FUNCTION

FUNCTION BadTimeArcUnit(BYVAL Unit AS LONG) AS LONG
SELECT CASE Unit
    CASE %UnitDegree
        FUNCTION = %FALSE
    CASE %UnitArcMin
        FUNCTION = %FALSE
    CASE %UnitArcSec
        FUNCTION = %FALSE
    CASE %UnitRadian
        FUNCTION = %FALSE
    CASE %UnitDay
        FUNCTION = %FALSE
    CASE %UnitHour
        FUNCTION = %FALSE
    CASE %UnitTMin
        FUNCTION = %FALSE
    CASE %UnitTSec
        FUNCTION = %FALSE
    CASE ELSE
        FUNCTION = %TRUE
END SELECT
END FUNCTION

FUNCTION BadTimeZone(BYVAL TimeZone AS DOUBLE) AS LONG

FUNCTION = (TimeZone >12. OR TimeZone < -13.)
END FUNCTION

FUNCTION BadStatistic(BYVAL Statistic AS LONG) AS LONG

SELECT CASE Statistic
    CASE %StatMaxPointMJD
        FUNCTION = %FALSE
    CASE %StatMaxPointSDY
        FUNCTION = %FALSE
    CASE %StatPrintAll
        FUNCTION = %FALSE
    CASE %StatMax
        FUNCTION = %FALSE
    CASE %StatAverage
        FUNCTION = %FALSE
    CASE ELSE
        FUNCTION = %TRUE

END SELECT


END FUNCTION






'END #INCLUDE "parmval.inc"      ' parameter validators


'#INCLUDE "astro.inc"        ' astronomical functions
 FUNCTION LH2LG(BYVAL LH AS EXT ) AS EXT
' Convert Heliocentric Longitude to Geocentric Longitude
' Meeus Astonomical Algorithms 2nd Ed. p166
FUNCTION = Reduce360(LH + 180)
END FUNCTION
 FUNCTION BH2Beta(BYVAL BH AS EXT  ) AS EXT
' Convert Heliocentric Latitude to Geocentric Latitude
' Meeus Astonomical Algorithms 2nd Ed. p166
FUNCTION = -BH
END FUNCTION
 FUNCTION LEB2Alpha(BYVAL xLambda  AS EXT , BYVAL xEpsilon  AS EXT , BYVAL xBeta AS EXT  ) AS EXT
' Right Ascension from Longitude, Obliquity and Latitude
' Coordinate Conversion
' Meeus Astonomical Algorithms 2nd Ed. p93
DIM  X AS EXT
DIM xLEB2Alpha AS EXT
X = dCos(xLambda)

xLEB2Alpha = dAtan((dSin(xLambda) * dCos(xEpsilon) - DTan(xBeta) * dSin(xEpsilon)) / X)
IF X < 0 THEN
    xLEB2Alpha = xLEB2Alpha + 180
END IF
FUNCTION = Reduce360(xLEB2Alpha)
END FUNCTION
FUNCTION LEB2Delta(BYVAL xLambda  AS EXT , BYVAL xEpsilon AS EXT  , BYVAL  xBeta AS EXT  ) AS EXT
' Geocentric Solar Declination
' Converstion of Units
' Meeus Astonomical Algorithms 2nd Ed. p93
FUNCTION = dASin(dSin(xBeta) * dCos(xEpsilon) + dCos(xBeta) * dSin(xEpsilon) * dSin(xLambda))
END FUNCTION

FUNCTION VSOP87Term(BYVAL VT  AS EXT, BYVAL A  AS EXT, BYVAL B  AS EXT, BYVAL C  AS EXT) AS EXT
FUNCTION = A * COS(B + C * VT)
END FUNCTION


 FUNCTION LGVSOP2FK5(BYVAL LG  AS EXT) AS EXT
' Convert Geocentric Longitude referred to dynamical mean equinox of date
' according to abridged VSOP87 TO FK5 reference system
' Meeus Astonomical Algorithms 2nd Ed. p166
FUNCTION = LG - DMS2Unit(0, 0, 0.09033,%UnitDegree)
END FUNCTION

FUNCTION EpsilonMeanIAU(BYVAL E AS EXT ) AS EXT
' Mean Obliquity of the Ecliptic
' Meeus Astonomical Algorithms 2nd Ed. p147

'KEYFILE V VALUES
'Professional Values
'Function = DMS2Unit(23, 26, 21.448,%UnitDegree) - _
'Suggested Freeware values
'Function = DMS2Unit(24, 17, 34.448,%UnitDegree) - _
FUNCTION = DMS2Unit(v(5), V(6), V(7),%UnitDegree) - _
    DMS2Unit(0, 0, 46.815 * E * 10,%UnitDegree) - _
    DMS2Unit(0, 0, 0.00059 * (E * 10) ^ 2,%UnitDegree) + _
    DMS2Unit(0, 0, 0.001813 * (E * 10) ^ 3,%UnitDegree) ' obliquity of the eliptic
'Professional Values
'Function = DMS2Unit(23, 26, 21.448,%UnitDegree) - _
'Suggested Freeware values
'Function = DMS2Unit(24, 17, 34.448,%UnitDegree) - _
END FUNCTION
 FUNCTION BetaVSOP2FK5(BYVAL xBeta AS EXT, BYVAL BV AS EXT, BYVAL LG AS EXT)    AS EXT
' Convert Geocentric Latitude referred to dynamical mean equinox of date
' according to abridged VSOP87 TO FK5 reference system
' Meeus Astonomical Algorithms 2nd Ed. p166
DIM  LambdaPrime AS EXT
LambdaPrime = LG - 1.397 * BV * 10 - 0.00031 * (BV * 10) ^ 2
FUNCTION = xBeta + DMS2Unit(0, 0, 0.03916,%UnitDegree) * (dCos(LambdaPrime) - dSin(LambdaPrime))
END FUNCTION
FUNCTION SunMeanLong(BYVAL SML  AS EXT) AS EXT
' Meeus Astonomical Algorithms 2nd Ed. p183
' According to VSOP87
FUNCTION = Reduce360(280.4664567 + 360007.6982779 * SML + 0.03032028 * SML ^ 2 + SML ^ 3 / 49931 - SML ^ 4 / 15300 - SML ^ 5 / 2000000)
END FUNCTION
 FUNCTION DIAU(BYVAL DI AS EXT ) AS EXT
' Mean elongation of the Moon from the Sun
' Meeus Astonomical Algorithms 2nd Ed. p144
FUNCTION = 297.85036 + 445267.11148 * DI * 10 - 0.0019142 * (DI * 10) ^ 2 + (DI * 10) ^ 3 / 189474
END FUNCTION
 FUNCTION MuIAU(BYVAL MI  AS EXT) AS EXT
' Mean anomaly of the Sun
' Meeus Astonomical Algorithms 2nd Ed. p144
FUNCTION = 357.52772 + 35999.05034 * MI * 10 - 0.0001603 * (MI * 10) ^ 2 - (MI * 10) ^ 3 / 300000
END FUNCTION
 FUNCTION MuPrimeIAU(BYVAL MPI AS EXT ) AS EXT
' Mean anomaly of the Moon
' Meeus Astonomical Algorithms 2nd Ed. p144
FUNCTION = 134.96298 + 477198.867398 * MPI * 10 + 0.0086972 * (MPI * 10) ^ 2 + (MPI * 10) ^ 3 / 56250
END FUNCTION
 FUNCTION FIAU(BYVAL FI AS EXT )  AS EXT
' Moon's argument of latitude
' Meeus Astonomical Algorithms 2nd Ed. p144
FUNCTION = 93.27191 + 483202.017538 * FI * 10 - 0.0036825 * (FI * 10) ^ 2 + (FI * 10) ^ 3 / 327270
END FUNCTION
 FUNCTION OmegaIAU(BYVAL OI AS EXT ) AS EXT
' Longitude of the ascENDing node of the moon's mean orbit on the ecliptic,
'   measured from the mean equinox of date
' Meeus Astonomical Algorithms 2nd Ed. p144
FUNCTION = Reduce360(125.04452 - 1934.136261 * OI * 10 + 0.0020708 * (OI * 10) ^ 2 + (OI * 10) ^ 3 / 450000)
END FUNCTION
 FUNCTION Omega3(BYVAL O  AS EXT)  AS EXT
' Longitude of the ascENDing node of the moon's mean orbit on the ecliptic,
'   measured from the mean equinox of date
' Meeus Astonomical Algorithms 2nd Ed. p144
FUNCTION = Reduce360(125.04452 - 1934.136261 * O * 10) '+ 0.0020708 * (O * 10) ^ 2 + (O * 10) ^ 3 / 450000)
END FUNCTION
 FUNCTION L3f(BYVAL L  AS EXT) AS EXT
' Mean Longitude of the Sun
' Meeus Astonomical Algorithms 2nd Ed. p144
FUNCTION = Reduce360(280.4665 + 36000.7698 * L * 10)
END FUNCTION
 FUNCTION LPrime3(BYVAL LP  AS EXT)  AS EXT
' Mean Longitude of the Moon
' Meeus Astonomical Algorithms 2nd Ed. p144
FUNCTION = Reduce360(218.3165 + 481267.8813 * LP * 10)
END FUNCTION
 FUNCTION DeltaPhi3(BYVAL DP AS EXT ) AS EXT
' Nutation in Longitude
' Meeus Astonomical Algorithms 2nd Ed. p144
' Corrected for removal of EXT negative sign before 2nd term
DIM DPOmega3 AS EXT
DPOmega3 = Omega3(DP)
FUNCTION = DMS2Unit(0, 0, -17.2,%UnitDegree) * dSin(DPOmega3) - _
    DMS2Unit(0, 0, 1.32,%UnitDegree) * _dSin(2 * L3f(DP)) - _
    DMS2Unit(0, 0, 0.23,%UnitDegree) * dSin(2 * LPrime3(DP)) + _
    DMS2Unit(0, 0, 0.21,%UnitDegree) * dSin(2 *  DPOmega3)
END FUNCTION
 FUNCTION DeltaEpsilon3(BYVAL DE AS EXT )  AS EXT
' Nutation in Obliquity
' Meeus Astonomical Algorithms 2nd Ed. p144
DIM DEOmega3 AS EXT
DEOmega3 = Omega3(DE)
FUNCTION = DMS2Unit(0, 0, 9.2,%UnitDegree) * dCos(DEOmega3) + DMS2Unit(0, 0, 0.57,%UnitDegree) * dCos(2 * L3f(DE)) + DMS2Unit(0, 0, 0.1,%UnitDegree) * dCos(2 * LPrime3(DE)) +DMS2Unit(0, 0, 0.09,%UnitDegree) * dCos(2 * DEOmega3)
END FUNCTION
 FUNCTION EccEarthOrbit(BYVAL EEO  AS EXT) AS EXT
' Meeus Astonomical Algorithms 2nd Ed. p163
FUNCTION = 0.016708634 - 0.000042037 * EEO * 10 - 0.0000001267 * (EEO * 10) ^ 2
END FUNCTION
 FUNCTION SunMeanAnom(BYVAL SMA AS EXT ) AS EXT
FUNCTION = Reduce360(357.52911 + 35999.05029 * SMA * 10 - 0.0001537 * (SMA * 10) ^ 2)
END FUNCTION
FUNCTION LHVSOP87a(BYVAL LV AS EXT)     AS EXT
'Earth's Heliocentric Ecliptical Longitude referred to dynamical mean equinox of date
' according to abridged VSOP87 in
' Meeus Astonomical Algorithms 2nd Ed. Ch32 and AppENDix III
DIM L0 AS EXT
DIM L1 AS EXT
DIM L2 AS EXT
DIM L3 AS EXT
DIM L4 AS EXT
DIM L5 AS EXT
DIM xLHVSOP87a AS EXT
L0 = 0

'KEYFILE V VALUES
'Model Lines with correct values
'L0 = L0 + VSOP87Term(LV, 175347046, 0, 0) '1
'L0 = L0 + VSOP87Term(LV, 3341656, 4.6692568, 6283.07585) '2
'Suggested Freeware values
'L0 = L0 + VSOP87Term(LV, 175407046, 0, 0) '1 in error
'L0 = L0 + VSOP87Term(LV, 3335656, 4.66, 6283.57585) '2 in error

L0 = L0 + VSOP87Term(LV, V(1), 0, 0) '1
L0 = L0 + VSOP87Term(LV, V(2), V(3), V(4)) '2
L0 = L0 + VSOP87Term(LV, 34894, 4.6261, 12566.1517) '3
L0 = L0 + VSOP87Term(LV, 3497, 2.7441, 5753.3849) '4
L0 = L0 + VSOP87Term(LV, 3418, 2.8289, 3.5231) '5

L0 = L0 + VSOP87Term(LV, 3136, 3.6277, 77713.7715) '6
L0 = L0 + VSOP87Term(LV, 2676, 4.4181, 7860.4194) '7
L0 = L0 + VSOP87Term(LV, 2343, 6.1352, 3930.2097) '8
L0 = L0 + VSOP87Term(LV, 1324, 0.7425, 11506.7698) '9
L0 = L0 + VSOP87Term(LV, 1273, 2.0371, 529.691) '10

L0 = L0 + VSOP87Term(LV, 1199, 1.1096, 1577.3435) '11
L0 = L0 + VSOP87Term(LV, 990, 5.233, 5884.927) '12
L0 = L0 + VSOP87Term(LV, 902, 2.045, 26.298) '13
L0 = L0 + VSOP87Term(LV, 857, 3.508, 398.149) '14
L0 = L0 + VSOP87Term(LV, 780, 1.179, 5223.694) '15

L0 = L0 + VSOP87Term(LV, 753, 2.533, 5507.553) '16
L0 = L0 + VSOP87Term(LV, 505, 4.583, 18849.228) '17
L0 = L0 + VSOP87Term(LV, 492, 4.205, 775.523) '18
L0 = L0 + VSOP87Term(LV, 357, 2.92, 0.067)  '19
L0 = L0 + VSOP87Term(LV, 317, 5.849, 11790.629) '20

L0 = L0 + VSOP87Term(LV, 284, 1.899, 796.298) '21
L0 = L0 + VSOP87Term(LV, 271, 0.315, 10977.079) '22
L0 = L0 + VSOP87Term(LV, 243, 0.345, 5486.778) '23
L0 = L0 + VSOP87Term(LV, 206, 4.806, 2544.314) '24
L0 = L0 + VSOP87Term(LV, 205, 1.869, 5573.143) '25

L0 = L0 + VSOP87Term(LV, 202, 2.458, 6069.777) '26
L0 = L0 + VSOP87Term(LV, 156, 0.833, 213.299) '27
L0 = L0 + VSOP87Term(LV, 132, 3.411, 2942.463) '28
L0 = L0 + VSOP87Term(LV, 126, 1.083, 20.775) '29
L0 = L0 + VSOP87Term(LV, 115, 0.645, 0.98)  '30

L0 = L0 + VSOP87Term(LV, 103, 0.636, 4694.003) '31
L0 = L0 + VSOP87Term(LV, 102, 0.976, 15720.839) '32
L0 = L0 + VSOP87Term(LV, 102, 4.267, 7.114) '33
L0 = L0 + VSOP87Term(LV, 99, 6.21, 2146.17)   '34
L0 = L0 + VSOP87Term(LV, 98, 0.68, 155.42)   '35

L0 = L0 + VSOP87Term(LV, 86, 5.98, 161000.69)   '36
L0 = L0 + VSOP87Term(LV, 85, 1.3, 6275.96)    '37
L0 = L0 + VSOP87Term(LV, 85, 3.67, 71430.7)    '38
L0 = L0 + VSOP87Term(LV, 80, 1.81, 17260.15)   '39
L0 = L0 + VSOP87Term(LV, 79, 3.04, 12036.46)   '40

L0 = L0 + VSOP87Term(LV, 75, 1.76, 5088.63)   '41
L0 = L0 + VSOP87Term(LV, 74, 3.5, 3154.69)    '42
L0 = L0 + VSOP87Term(LV, 74, 4.68, 801.82)   '43
L0 = L0 + VSOP87Term(LV, 70, 0.83, 9437.76)   '44
L0 = L0 + VSOP87Term(LV, 62, 3.98, 8827.39)   '45

L0 = L0 + VSOP87Term(LV, 61, 1.82, 7084.9)    '46
L0 = L0 + VSOP87Term(LV, 57, 2.78, 6286.6)    '47
L0 = L0 + VSOP87Term(LV, 56, 4.39, 14143.5)    '48
L0 = L0 + VSOP87Term(LV, 56, 3.47, 6279.55)   '49
L0 = L0 + VSOP87Term(LV, 52, 0.19, 12139.55)   '50

' Tried here not accurate enough
L0 = L0 + VSOP87Term(LV, 52, 1.33, 1748.02)   '51
L0 = L0 + VSOP87Term(LV, 51, 0.28, 5856.48)   '52
L0 = L0 + VSOP87Term(LV, 49, 0.49, 1194.45)   '53
L0 = L0 + VSOP87Term(LV, 41, 5.37, 8429.24)   '54
L0 = L0 + VSOP87Term(LV, 41, 2.4, 19651.05)    '55

L0 = L0 + VSOP87Term(LV, 39, 6.17, 10447.39)   '56
L0 = L0 + VSOP87Term(LV, 37, 6.04, 10213.29)   '57
L0 = L0 + VSOP87Term(LV, 37, 2.57, 1059.38)   '58
L0 = L0 + VSOP87Term(LV, 36, 1.71, 2352.87)   '59
L0 = L0 + VSOP87Term(LV, 36, 1.78, 6812.77)   '60

' Tried here not accurate enough
L0 = L0 + VSOP87Term(LV, 33, 0.59, 17789.85)   '61
L0 = L0 + VSOP87Term(LV, 30, 0.44, 83996.85)   '62
L0 = L0 + VSOP87Term(LV, 30, 2.74, 1349.87)   '63
L0 = L0 + VSOP87Term(LV, 25, 3.16, 4690.48)   '64

L1 = 0
L1 = L1 + VSOP87Term(LV, 628331966747#, 0, 0) '1


L1 = L1 + VSOP87Term(LV, 206059, 2.678235, 6283.07585)  '2
L1 = L1 + VSOP87Term(LV, 4303, 2.6351, 12566.1517) '3

L1 = L1 + VSOP87Term(LV, 425, 1.59, 3.523) '4
L1 = L1 + VSOP87Term(LV, 119, 5.796, 26.298) '5

' Tried here - too accurate
L1 = L1 + VSOP87Term(LV, 109, 2.966, 1577.334) '6
L1 = L1 + VSOP87Term(LV, 93, 2.59, 18849.23)   '7
L1 = L1 + VSOP87Term(LV, 72, 1.14, 529.69)   '8
L1 = L1 + VSOP87Term(LV, 68, 1.87, 398.15)   '9
L1 = L1 + VSOP87Term(LV, 67, 4.41, 5507.55)   '10

' Tried here too accurate
L1 = L1 + VSOP87Term(LV, 59, 2.89, 5223.69)   '11
L1 = L1 + VSOP87Term(LV, 56, 2.17, 155.42)   '12
L1 = L1 + VSOP87Term(LV, 45, 0.4, 796.3)     '13
L1 = L1 + VSOP87Term(LV, 36, 0.47, 775.52)   '14
L1 = L1 + VSOP87Term(LV, 29, 2.65, 7.11)   '15

L1 = L1 + VSOP87Term(LV, 21, 5.34, 0.98)   '16
L1 = L1 + VSOP87Term(LV, 19, 1.85, 5486.78)   '17
L1 = L1 + VSOP87Term(LV, 19, 4.97, 213.3)    '18
L1 = L1 + VSOP87Term(LV, 17, 2.99, 6275.96)   '19
L1 = L1 + VSOP87Term(LV, 16, 0.03, 2544.31)   '20

L1 = L1 + VSOP87Term(LV, 16, 1.43, 2146.17)   '21
L1 = L1 + VSOP87Term(LV, 15, 1.21, 10977.08)   '22
L1 = L1 + VSOP87Term(LV, 12, 2.83, 1748.02)   '23
L1 = L1 + VSOP87Term(LV, 12, 3.26, 5088.63)   '24
L1 = L1 + VSOP87Term(LV, 12, 5.27, 1194.45)   '25

L1 = L1 + VSOP87Term(LV, 12, 2.08, 4694#)     '26
L1 = L1 + VSOP87Term(LV, 11, 0.77, 553.57)   '27
L1 = L1 + VSOP87Term(LV, 10, 1.3, 6286.6)     '28
L1 = L1 + VSOP87Term(LV, 10, 4.24, 1349.87)   '29
L1 = L1 + VSOP87Term(LV, 9, 2.7, 242.73)    '30

L1 = L1 + VSOP87Term(LV, 9, 5.64, 951.72)   '31
L1 = L1 + VSOP87Term(LV, 8, 5.3, 2352.87)    '32
L1 = L1 + VSOP87Term(LV, 6, 2.65, 9437.76)   '33
L1 = L1 + VSOP87Term(LV, 6, 4.67, 4690.48)   '34

L2 = 0
L2 = L2 + VSOP87Term(LV, 52919, 0, 0) '1
L2 = L2 + VSOP87Term(LV, 8720, 1.0721, 6283.0758) '2
L2 = L2 + VSOP87Term(LV, 309, 0.867, 12566.152) '3
L2 = L2 + VSOP87Term(LV, 27, 0.05, 3.52) '4
L2 = L2 + VSOP87Term(LV, 16, 5.19, 26.3) '5

L2 = L2 + VSOP87Term(LV, 16, 3.68, 155.42)   '6
L2 = L2 + VSOP87Term(LV, 10, 0.76, 18849.23)   '7
L2 = L2 + VSOP87Term(LV, 9, 2.06, 77713.77)   '8
L2 = L2 + VSOP87Term(LV, 7, 0.83, 775.52)   '9
L2 = L2 + VSOP87Term(LV, 5, 4.66, 1577.34)   '10

L2 = L2 + VSOP87Term(LV, 4, 1.03, 7.11)   '11
L2 = L2 + VSOP87Term(LV, 4, 3.44, 5573.14)   '12
L2 = L2 + VSOP87Term(LV, 3, 5.14, 796.3)    '13
L2 = L2 + VSOP87Term(LV, 3, 6.05, 5507.55)   '14
L2 = L2 + VSOP87Term(LV, 3, 1.19, 242.73)   '15

L2 = L2 + VSOP87Term(LV, 3, 6.12, 529.69)   '16
L2 = L2 + VSOP87Term(LV, 3, 0.31, 398.15)   '17
L2 = L2 + VSOP87Term(LV, 3, 2.28, 553.57)   '18
L2 = L2 + VSOP87Term(LV, 2, 4.38, 5223.69)   '19
L2 = L2 + VSOP87Term(LV, 2, 3.75, 0.98)   '20

L3 = 0
L3 = L3 + VSOP87Term(LV, 289, 5.844, 6283.076) '1
L3 = L3 + VSOP87Term(LV, 35, 0, 0) '2
L3 = L3 + VSOP87Term(LV, 17, 5.49, 12566.15) '3
L3 = L3 + VSOP87Term(LV, 3, 5.2, 155.42) '4
L3 = L3 + VSOP87Term(LV, 1, 4.72, 3.52) '5

L3 = L3 + VSOP87Term(LV, 1, 5.3, 18849.23) '6
L3 = L3 + VSOP87Term(LV, 1, 5.97, 242.73) '7 complete

L4 = 0
L4 = L4 + VSOP87Term(LV, 144, 3.142, 0) '1
L4 = L4 + VSOP87Term(LV, 8, 4.13, 6283.08) '2
L4 = L4 + VSOP87Term(LV, 1, 3.84, 12566.15) '3 complete
L5 = VSOP87Term(LV, 1, 3.14, 0) '1 complete


' In radians

xLHVSOP87a = (L0 + L1 * LV + L2 * LV ^ 2 + L3 * LV ^ 3 + L4 * LV ^ 4 + L5 * LV ^ 5) / 10 ^ 8

' In degrees between [0,360)
' && round to 1/1000000 of a degree
FUNCTION = Reduce360(xLHVSOP87a * 57.2957795130823)

END FUNCTION
FUNCTION RVSOP87a(BYVAL RV AS EXT)  AS EXT
'Earth's orbital Radius vector in AU
' according to abridged VSOP87 in
' Meeus Astonomical Algorithms 2nd Ed. Ch32 and AppENDix III
DIM R0 AS EXT
DIM R1 AS EXT
DIM R2 AS EXT
DIM R3 AS EXT
DIM R4 AS EXT

R0 = 0
R0 = R0 + VSOP87Term(RV, 100013989, 0#, 0#)            '1
R0 = R0 + VSOP87Term(RV, 1670700, 3.0984635, 6283.07585) '2
R0 = R0 + VSOP87Term(RV, 13956, 3.05525, 12566.1517)   '3
R0 = R0 + VSOP87Term(RV, 3084, 5.1985, 77713.7715)    '4
R0 = R0 + VSOP87Term(RV, 1628, 1.1739, 5753.3849)    '5


R0 = R0 + VSOP87Term(RV, 1576, 2.8469, 7860.4194)    '6
R0 = R0 + VSOP87Term(RV, 925, 5.453, 11506.77)       '7
R0 = R0 + VSOP87Term(RV, 542, 4.564, 3930.21)       '8
R0 = R0 + VSOP87Term(RV, 472, 3.661, 5884.927)      '9
R0 = R0 + VSOP87Term(RV, 346, 0.964, 5507.553)      '10

R0 = R0 + VSOP87Term(RV, 329, 5.9, 5223.694)        '11
R0 = R0 + VSOP87Term(RV, 307, 0.299, 5573.143)      '12
R0 = R0 + VSOP87Term(RV, 243, 4.273, 11790.629)      '13
R0 = R0 + VSOP87Term(RV, 212, 5.847, 1577.344)      '14
R0 = R0 + VSOP87Term(RV, 186, 5.022, 10977.079)      '15

R0 = R0 + VSOP87Term(RV, 175, 3.012, 18849.228)      '16
R0 = R0 + VSOP87Term(RV, 110, 5.055, 5486.778)      '17
R0 = R0 + VSOP87Term(RV, 98, 0.89, 6069.78)        '18
R0 = R0 + VSOP87Term(RV, 86, 5.69, 15720.84)        '19
R0 = R0 + VSOP87Term(RV, 86, 1.27, 161000.69)        '20

R0 = R0 + VSOP87Term(RV, 65, 0.27, 17260.15)        '21
R0 = R0 + VSOP87Term(RV, 63, 0.92, 529.69)        '22
R0 = R0 + VSOP87Term(RV, 57, 2.01, 83996.85)        '23
R0 = R0 + VSOP87Term(RV, 56, 5.24, 71430.7)         '24
R0 = R0 + VSOP87Term(RV, 49, 3.25, 2544.31)        '25

R0 = R0 + VSOP87Term(RV, 47, 2.58, 775.52)        '26
R0 = R0 + VSOP87Term(RV, 45, 5.54, 9437.76)        '27
R0 = R0 + VSOP87Term(RV, 43, 6.01, 6275.96)        '28
R0 = R0 + VSOP87Term(RV, 39, 5.36, 4694#)          '29
R0 = R0 + VSOP87Term(RV, 38, 2.39, 8827.39)        '30

R0 = R0 + VSOP87Term(RV, 37, 0.83, 19651.05)        '31
R0 = R0 + VSOP87Term(RV, 37, 4.9, 12139.55)         '32
R0 = R0 + VSOP87Term(RV, 36, 1.67, 12036.46)        '33
R0 = R0 + VSOP87Term(RV, 35, 1.84, 2942.46)        '34
R0 = R0 + VSOP87Term(RV, 33, 0.24, 7084.9)         '35

R0 = R0 + VSOP87Term(RV, 32, 0.18, 5088.63)        '36
R0 = R0 + VSOP87Term(RV, 32, 1.78, 398.15)        '37
R0 = R0 + VSOP87Term(RV, 28, 1.21, 6286.6)         '38
R0 = R0 + VSOP87Term(RV, 28, 1.9, 6279.55)         '39
R0 = R0 + VSOP87Term(RV, 26, 4.59, 10447.39)        '40

R1 = 0
R1 = R1 + VSOP87Term(RV, 103019, 1.10749, 6283.07585)  '1
R1 = R1 + VSOP87Term(RV, 1721, 1.0644, 12566.1517)    '2
R1 = R1 + VSOP87Term(RV, 702, 3.142, 0#)         '3
R1 = R1 + VSOP87Term(RV, 32, 1.02, 18849.23)        '4
R1 = R1 + VSOP87Term(RV, 31, 2.84, 5507.55)        '5

R1 = R1 + VSOP87Term(RV, 25, 1.32, 5223.69)        '6
R1 = R1 + VSOP87Term(RV, 18, 1.42, 1577.34)        '7
R1 = R1 + VSOP87Term(RV, 10, 5.91, 10977.08)        '8
R1 = R1 + VSOP87Term(RV, 9, 1.42, 6275.96)        '9
R1 = R1 + VSOP87Term(RV, 9, 0.27, 5486.78)        '10

R2 = 0
R2 = R2 + VSOP87Term(RV, 4359, 5.7846, 6283.0758)    '1
R2 = R2 + VSOP87Term(RV, 124, 5.579, 12566.152)      '2
R2 = R2 + VSOP87Term(RV, 12, 3.14, 0#)          '3
R2 = R2 + VSOP87Term(RV, 9, 3.63, 77713.77)        '4
R2 = R2 + VSOP87Term(RV, 6, 1.87, 5573.14)        '5

R2 = R2 + VSOP87Term(RV, 3, 5.47, 18849.23)        '6

R3 = 0
R3 = R3 + VSOP87Term(RV, 145, 4.273, 6283.076)        '1
R3 = R3 + VSOP87Term(RV, 7, 3.92, 12566.15)      '2

R4 = 0
R4 = R4 + VSOP87Term(RV, 4, 2.56, 6283.08)    '1

FUNCTION = (R0 + R1 * RV + R2 * RV ^ 2 + R3 * RV ^ 3 + R4 * RV ^ 4) / 10 ^ 8
END FUNCTION
FUNCTION BHVSOP87a(BYVAL BV AS EXT) AS EXT
'Earth's Heliocentric Ecliptical Latitude referred to dynamical mean equinox of date
' according to abridged VSOP87 in
' Meeus Astonomical Algorithms 2nd Ed. Ch32 and AppENDix III
DIM B0 AS EXT
DIM B1 AS EXT
DIM xBHVSOP87a AS EXT
B0 = 0
B0 = B0 + VSOP87Term(BV, 280, 3.199, 84334.662) '1

B0 = B0 + VSOP87Term(BV, 102, 5.422, 5507.553) '2
B0 = B0 + VSOP87Term(BV, 80, 3.88, 5223.69) '3
B0 = B0 + VSOP87Term(BV, 44, 3.7, 2352.87) '4
B0 = B0 + VSOP87Term(BV, 32, 4, 1577.34) '5

B1 = 0
B1 = B1 + VSOP87Term(BV, 9, 3.9, 5507.55) '1
B1 = B1 + VSOP87Term(BV, 6, 1.73, 5223.69) '2

' In radians

xBHVSOP87a = (B0 + B1 * BV) / 10 ^ 8

' In degrees between [0,360)
' && round to 1/1000000 of a degree
FUNCTION = xBHVSOP87a * 57.2957795130823
END FUNCTION

FUNCTION DeltaPhiIAU(BYVAL DPI AS EXT)  AS EXT
' Nutation in Longitude IAU 1984
' Meeus Astonomical Algorithms 2nd Ed. p144-6
DIM d AS EXT
DIM Mu AS EXT
DIM MuPrime AS EXT
DIM F AS EXT
DIM Omega AS EXT
DIM xDeltaPhiIAU AS EXT
d = DIAU(DPI)
Mu = MuIAU(DPI)
MuPrime = MuPrimeIAU(DPI)
F = FIAU(DPI)
Omega = OmegaIAU(DPI)


xDeltaPhiIAU = 0

' Unit for DeltaPhiIAU during these periodic terms is 10000 arcsecond

xDeltaPhiIAU = xDeltaPhiIAU + (-171996 - 174.2 * DPI * 10) * dSin(Omega) '1
xDeltaPhiIAU = xDeltaPhiIAU + (-13187 - 1.6 * DPI * 10) * dSin(-2 * d + 2 * F + 2 * Omega) '2
xDeltaPhiIAU = xDeltaPhiIAU + (-2274 - 0.2 * DPI * 10) * dSin(2 * F + 2 * Omega) ' 3
xDeltaPhiIAU = xDeltaPhiIAU + (2062 + 0.2 * DPI * 10) * dSin(2 * Omega) ' 4
xDeltaPhiIAU = xDeltaPhiIAU + (1426 - 3.4 * DPI * 10) * dSin(Mu) ' 5

xDeltaPhiIAU = xDeltaPhiIAU + (712 + 0.1 * DPI * 10) * dSin(MuPrime) ' 6
xDeltaPhiIAU = xDeltaPhiIAU + (-517 + 1.2 * DPI * 10) * dSin(-2 * d + Mu + 2 * F + 2 * Omega) ' 7
xDeltaPhiIAU = xDeltaPhiIAU + (-386 - 0.4 * DPI * 10) * dSin(2 * F + Omega) ' 8
xDeltaPhiIAU = xDeltaPhiIAU - 301 * dSin(MuPrime + 2 * F + 2 * Omega) ' 9
xDeltaPhiIAU = xDeltaPhiIAU + (217 - 0.5 * DPI * 10) * dSin(-2 * d - Mu + 2 * F + 2 * Omega) ' 10

xDeltaPhiIAU = xDeltaPhiIAU - 158 * dSin(-2 * d + MuPrime) ' 11
xDeltaPhiIAU = xDeltaPhiIAU + (129 + 0.1 * DPI * 10) * dSin(-2 * d + 2 * F + Omega) ' 12
xDeltaPhiIAU = xDeltaPhiIAU + 123 * dSin(-1 * MuPrime + 2 * F + 2 * Omega) ' 13
xDeltaPhiIAU = xDeltaPhiIAU + 63 * dSin(2 * d)  ' 14
xDeltaPhiIAU = xDeltaPhiIAU + (63 + 0.1 * DPI * 10) * dSin(MuPrime + Omega) ' 15

xDeltaPhiIAU = xDeltaPhiIAU - 59 * dSin(2 * d - MuPrime + 2 * F + 2 * Omega) ' 16
xDeltaPhiIAU = xDeltaPhiIAU + (-58 - 0.1 * DPI * 10) * dSin(-1 * MuPrime + Omega) ' 17
xDeltaPhiIAU = xDeltaPhiIAU - 51 * dSin(MuPrime + 2 * F + Omega) ' 18
xDeltaPhiIAU = xDeltaPhiIAU + 48 * dSin(-2 * d + 2 * MuPrime) ' 19
xDeltaPhiIAU = xDeltaPhiIAU + 46 * dSin(-2 * MuPrime + 2 * F + Omega) ' 20

xDeltaPhiIAU = xDeltaPhiIAU - 38 * dSin(2 * d + 2 * F + 2 * Omega) ' 21
xDeltaPhiIAU = xDeltaPhiIAU - 31 * dSin(2 * MuPrime + 2 * F + 2 * Omega) ' 22
xDeltaPhiIAU = xDeltaPhiIAU + 29 * dSin(2 * MuPrime) ' 23
xDeltaPhiIAU = xDeltaPhiIAU + 29 * dSin(-2 * d + MuPrime + 2 * F + 2 * Omega) ' 24
xDeltaPhiIAU = xDeltaPhiIAU + 26 * dSin(2 * F) ' 25

xDeltaPhiIAU = xDeltaPhiIAU - 22 * dSin(-2 * d + 2 * F)  ' 26
xDeltaPhiIAU = xDeltaPhiIAU + 21 * dSin(-1 * MuPrime + 2 * F + Omega) ' 27
xDeltaPhiIAU = xDeltaPhiIAU + (17 - 0.1 * DPI * 10) * dSin(2 * Mu) ' 28
xDeltaPhiIAU = xDeltaPhiIAU + 16 * dSin(2 * d - MuPrime + Omega)   ' 29
xDeltaPhiIAU = xDeltaPhiIAU + (-16 + 0.1 * DPI * 10) * dSin(-2 * d + 2 * Mu + 2 * F + 2 * Omega) ' 30

xDeltaPhiIAU = xDeltaPhiIAU - 15 * dSin(Mu + Omega)  ' 31
xDeltaPhiIAU = xDeltaPhiIAU - 13 * dSin(-2 * d + MuPrime + Omega) ' 32
xDeltaPhiIAU = xDeltaPhiIAU - 12 * dSin(-1 * Mu + Omega) ' 33
xDeltaPhiIAU = xDeltaPhiIAU + 11 * dSin(2 * MuPrime - 2 * F) ' 34
xDeltaPhiIAU = xDeltaPhiIAU - 10 * dSin(2 * d - MuPrime + 2 * F + Omega) ' 35

xDeltaPhiIAU = xDeltaPhiIAU - 8 * dSin(2 * d + MuPrime + 2 * F + 2 * Omega) ' 36
xDeltaPhiIAU = xDeltaPhiIAU + 7 * dSin(Mu + 2 * F + 2 * Omega) ' 37
xDeltaPhiIAU = xDeltaPhiIAU - 7 * dSin(-2 * d + Mu + MuPrime) ' 38
xDeltaPhiIAU = xDeltaPhiIAU - 7 * dSin(-1 * Mu + 2 * F + 2 * Omega) ' 39
xDeltaPhiIAU = xDeltaPhiIAU - 7 * dSin(2 * d + 2 * F + Omega) ' 40

xDeltaPhiIAU = xDeltaPhiIAU + 6 * dSin(2 * d + MuPrime)  ' 41
xDeltaPhiIAU = xDeltaPhiIAU + 6 * dSin(-2 * d + 2 * MuPrime + 2 * F + 2 * Omega) ' 42
xDeltaPhiIAU = xDeltaPhiIAU + 6 * dSin(-2 * d + MuPrime + 2 * F + Omega) ' 43
xDeltaPhiIAU = xDeltaPhiIAU - 6 * dSin(2 * d - 2 * MuPrime + Omega)  ' 44
xDeltaPhiIAU = xDeltaPhiIAU - 6 * dSin(2 * d + Omega)   ' 45

xDeltaPhiIAU = xDeltaPhiIAU + 5 * dSin(-1 * Mu + MuPrime) ' 46
xDeltaPhiIAU = xDeltaPhiIAU - 5 * dSin(-2 * d - Mu + 2 * F + Omega) ' 47
xDeltaPhiIAU = xDeltaPhiIAU - 5 * dSin(-2 * d + Omega)   ' 48
xDeltaPhiIAU = xDeltaPhiIAU - 5 * dSin(2 * MuPrime + 2 * F + Omega) ' 49
xDeltaPhiIAU = xDeltaPhiIAU + 4 * dSin(-2 * d + 2 * MuPrime + Omega)  ' 50

xDeltaPhiIAU = xDeltaPhiIAU + 4 * dSin(-2 * d + Mu + 2 * F + Omega)  ' 51
xDeltaPhiIAU = xDeltaPhiIAU + 4 * dSin(MuPrime - 2 * F) ' 52
xDeltaPhiIAU = xDeltaPhiIAU - 4 * dSin(-1 * d + MuPrime)   ' 53
xDeltaPhiIAU = xDeltaPhiIAU - 4 * dSin(-2 * d + Mu)   ' 54
xDeltaPhiIAU = xDeltaPhiIAU - 4 * dSin(d)   ' 55

xDeltaPhiIAU = xDeltaPhiIAU + 3 * dSin(MuPrime + 2 * F) ' 56
xDeltaPhiIAU = xDeltaPhiIAU - 3 * dSin(-2 * MuPrime + 2 * F + 2 * Omega) ' 57
xDeltaPhiIAU = xDeltaPhiIAU - 3 * dSin(-1 * d - Mu + MuPrime) ' 58
xDeltaPhiIAU = xDeltaPhiIAU - 3 * dSin(Mu + MuPrime) ' 59
xDeltaPhiIAU = xDeltaPhiIAU - 3 * dSin(-1 * Mu + MuPrime + 2 * F + 2 * Omega) ' 60

xDeltaPhiIAU = xDeltaPhiIAU - 3 * dSin(2 * d - Mu - MuPrime + 2 * F + 2 * Omega) ' 61
xDeltaPhiIAU = xDeltaPhiIAU - 3 * dSin(3 * MuPrime + 2 * F + 2 * Omega) ' 62
xDeltaPhiIAU = xDeltaPhiIAU - 3 * dSin(2 * d - Mu + 2 * F + 2 * Omega) ' 63
' Convert to arc seconds and round to 1/1000 arc second accuracy.
' minimum error is +/- 3/10000 arc seconds therefore round to approx. nearest 1/1000 of arcsecond
' after Converting to degrees
FUNCTION = DMS2Unit(0, 0, xDeltaPhiIAU / 10000,%UnitDegree)
END FUNCTION
FUNCTION DeltaEpsilonIAU(BYVAL DEI  AS EXT)  AS EXT
' Nutation in Obliquity IAU 1984
' Meeus Astonomical Algorithms 2nd Ed. p144-6

DIM d AS EXT
DIM Mu AS EXT
DIM MuPrime AS EXT
DIM F AS EXT
DIM Omega AS EXT
DIM xDeltaEpsilonIAU AS EXT

d = DIAU(DEI)
Mu = MuIAU(DEI)
MuPrime = MuPrimeIAU(DEI)
F = FIAU(DEI)
Omega = OmegaIAU(DEI)


xDeltaEpsilonIAU = 0

' Unit for xDeltaEpsilonIAU during these periodic terms is 10000 arcsecond

xDeltaEpsilonIAU = xDeltaEpsilonIAU + (92025 + 8.9 * DEI * 10) * dCos(Omega) '1
xDeltaEpsilonIAU = xDeltaEpsilonIAU + (5736 - 3.1 * DEI * 10) * dCos(-2 * d + 2 * F + 2 * Omega) '2
xDeltaEpsilonIAU = xDeltaEpsilonIAU + (977 - 0.5 * DEI * 10) * dCos(2 * F + 2 * Omega) ' 3
xDeltaEpsilonIAU = xDeltaEpsilonIAU + (-895 + 0.5 * DEI * 10) * dCos(2 * Omega) ' 4
xDeltaEpsilonIAU = xDeltaEpsilonIAU + (54 - 0.1 * DEI * 10) * dCos(Mu) ' 5

xDeltaEpsilonIAU = xDeltaEpsilonIAU - 7 * dCos(MuPrime) ' 6
xDeltaEpsilonIAU = xDeltaEpsilonIAU + (224 - 0.6 * DEI * 10) * dCos(-2 * d + Mu + 2 * F + 2 * Omega) ' 7
xDeltaEpsilonIAU = xDeltaEpsilonIAU + 200 * dCos(2 * F + Omega) ' 8
xDeltaEpsilonIAU = xDeltaEpsilonIAU + (129 - 0.1 * DEI * 10) * dCos(MuPrime + 2 * F + 2 * Omega) ' 9
xDeltaEpsilonIAU = xDeltaEpsilonIAU + (-95 + 0.3 * DEI * 10) * dCos(-2 * d - Mu + 2 * F + 2 * Omega) ' 10

' 11
xDeltaEpsilonIAU = xDeltaEpsilonIAU - 70 * dCos(-2 * d + 2 * F + Omega) ' 12
xDeltaEpsilonIAU = xDeltaEpsilonIAU - 53 * dCos(-1 * MuPrime + 2 * F + 2 * Omega) ' 13
' 14
xDeltaEpsilonIAU = xDeltaEpsilonIAU + (-33) * dCos(MuPrime + Omega) ' 15

xDeltaEpsilonIAU = xDeltaEpsilonIAU + 26 * dCos(2 * d - MuPrime + 2 * F + 2 * Omega) ' 16
xDeltaEpsilonIAU = xDeltaEpsilonIAU + 32 * dCos(-1 * MuPrime + Omega) ' 17
xDeltaEpsilonIAU = xDeltaEpsilonIAU + 27 * dCos(MuPrime + 2 * F + Omega) ' 18
' 19
xDeltaEpsilonIAU = xDeltaEpsilonIAU - 24 * dCos(-2 * MuPrime + 2 * F + Omega) ' 20

xDeltaEpsilonIAU = xDeltaEpsilonIAU + 16 * dCos(2 * d + 2 * F + 2 * Omega) ' 21
xDeltaEpsilonIAU = xDeltaEpsilonIAU + 13 * dCos(2 * MuPrime + 2 * F + 2 * Omega) ' 22
' 23
xDeltaEpsilonIAU = xDeltaEpsilonIAU - 12 * dCos(-2 * d + MuPrime + 2 * F + 2 * Omega) ' 24
' 25

' 26
xDeltaEpsilonIAU = xDeltaEpsilonIAU - 10 * dCos(-1 * MuPrime + 2 * F + Omega) ' 27
' 28
xDeltaEpsilonIAU = xDeltaEpsilonIAU - 8 * dCos(2 * d - MuPrime + Omega)   ' 29
xDeltaEpsilonIAU = xDeltaEpsilonIAU + 7 * dCos(-2 * d + 2 * Mu + 2 * F + 2 * Omega) ' 30

xDeltaEpsilonIAU = xDeltaEpsilonIAU + 9 * dCos(Mu + Omega) ' 31
xDeltaEpsilonIAU = xDeltaEpsilonIAU + 7 * dCos(-2 * d + MuPrime + Omega) ' 32
xDeltaEpsilonIAU = xDeltaEpsilonIAU + 6 * dCos(-1 * Mu + Omega) ' 33
' 34
xDeltaEpsilonIAU = xDeltaEpsilonIAU + 5 * dCos(2 * d - MuPrime + 2 * F + Omega) ' 35

xDeltaEpsilonIAU = xDeltaEpsilonIAU + 3 * dCos(2 * d + MuPrime + 2 * F + 2 * Omega) ' 36
xDeltaEpsilonIAU = xDeltaEpsilonIAU - 3 * dCos(Mu + 2 * F + 2 * Omega) ' 37
' 38
xDeltaEpsilonIAU = xDeltaEpsilonIAU + 3 * dCos(-1 * Mu + 2 * F + 2 * Omega) ' 39
xDeltaEpsilonIAU = xDeltaEpsilonIAU + 3 * dCos(2 * d + 2 * F + Omega) ' 40

' 41
xDeltaEpsilonIAU = xDeltaEpsilonIAU - 3 * dCos(-2 * d + 2 * MuPrime + 2 * F + 2 * Omega) ' 42
xDeltaEpsilonIAU = xDeltaEpsilonIAU - 3 * dCos(-2 * d + MuPrime + 2 * F + Omega) ' 43
xDeltaEpsilonIAU = xDeltaEpsilonIAU + 3 * dCos(2 * d - 2 * MuPrime + Omega) ' 44
xDeltaEpsilonIAU = xDeltaEpsilonIAU + 3 * dCos(2 * d + Omega)  ' 45

' 46
xDeltaEpsilonIAU = xDeltaEpsilonIAU + 3 * dCos(-2 * d - Mu + 2 * F + Omega) ' 47
xDeltaEpsilonIAU = xDeltaEpsilonIAU + 3 * dCos(-2 * d + Omega)  ' 48
xDeltaEpsilonIAU = xDeltaEpsilonIAU + 3 * dCos(2 * MuPrime + 2 * F + Omega) ' 49

' Convert to arc seconds and round to 1/1000 arc second accuracy.
' minimum error is +/-  3/10000 arc seconds therefore round to approx. nearest 1/1000 of arcsecond
' after Converting to degrees
FUNCTION = DMS2Unit(0, 0, xDeltaEpsilonIAU / 10000,%UnitDegree)

END FUNCTION

FUNCTION Epsilon(BYVAL EE AS EXT ) AS EXT
' Apparent Obliquity of the Ecliptic
' Meeus Astonomical Algorithms 2nd Ed. p147
'Epsilon = EpsilonMeanLaskar(EE*10) + DeltaEpsilonIAU(EE*10) ' Laskar is more accurate for mean
FUNCTION = EpsilonMeanIAU(EE) + DeltaEpsilonIAU(EE)  ' but IAU FUNCTIONs work together

END FUNCTION
 FUNCTION Beta(BYVAL BF AS EXT) AS EXT
'Sun's Geocentric Ecliptical Latitude referred to FK5
' Meeus Astonomical Algorithms 2nd Ed. p166
FUNCTION = BetaVSOP2FK5(BH2Beta(BHVSOP87a(BF)), BF, LH2LG(LHVSOP87a(BF)))
END FUNCTION
 FUNCTION LTrue(BYVAL LT AS EXT)AS EXT
'Sun's True Geocentric Ecliptical Longitude referred to FK5
' Meeus Astonomical Algorithms 2nd Ed. p166
FUNCTION = LGVSOP2FK5(LH2LG(LHVSOP87a(LT)))
END FUNCTION


FUNCTION Lambda(BYVAL LA AS EXT ) AS EXT
'Sun's Apparent Geocentric Ecliptical Longitude referred to FK5
' after accounting for nutation and aberration
' Meeus Astonomical Algorithms 2nd Ed. p166
DIM xLambda AS EXT
xLambda = LTrue(LA) + DeltaPhiIAU(LA)  ' nutation adjustment
xLambda = xLambda - DMS2Unit(0, 0, 20.4898,%UnitDegree) / RVSOP87a(LA) ' aberration adjustment
FUNCTION = xLambda
END FUNCTION


FUNCTION Alpha(BYVAL AA AS EXT )   AS EXT
' Geocentric Solar Right Ascension
' Meeus Astonomical Algorithms 2nd Ed. p166-68
FUNCTION = LEB2Alpha(Lambda(AA), Epsilon(AA), Beta(AA))
END FUNCTION


FUNCTION Delta(BYVAL D AS EXT ) AS EXT
' Solar Declination
' from VSOP87
' Meeus Astonomical Algorithms 2nd Ed. p166-68
FUNCTION = LEB2Delta(Lambda(D), Epsilon(D), Beta(D))
END FUNCTION

FUNCTION EOT(BYVAL E AS EXT ) AS EXT
' Equation of Time
' High Accuracy
' Meeus p 183
FUNCTION = Reduce180Signed(SunMeanLong(E) - 0.0057183 - LEB2Alpha(Lambda(E), Epsilon(E), Beta(E)) + DeltaPhiIAU(E) * dCos(Epsilon(E))) '* 4  - 0.20 / 60 '0.2 sec adjustment to track MICA

END FUNCTION


FUNCTION Sun (BYVAL B AS EXT, BYVAL SIndex AS LONG) AS EXT

'STATIC EXT V[]
'Remove Static element
'Now the function does extra calculation for each call.  Could be streamlined.
'or abandoned and use the individual functions instead.
'
'EXT V[]

DIM CEOT AS EXT
DIM CAlpha         AS EXT
DIM CLtrue                   AS EXT
DIM CDeltaPhi                          AS EXT
DIM CDeltaEpsilon                                AS EXT
DIM CSunMeanLong                                           AS EXT
DIM CDelta AS EXT
DIM CEpsilonMean     AS EXT
DIM CRvsop87a                  AS EXT
DIM CBhvsop87a                           AS EXT
DIM CLhvsop87a                                     AS EXT
DIM CLG                                                      AS EXT
DIM CBeta AS EXT
DIM CLambda         AS EXT
DIM CEpsilon                  AS EXT
DIM Cwidth AS LONG


'''VSize = 18

' Remove statid elements
'IFZ V[] THEN
'''     DIM V[VSize] ' Dimension on first entry
'''     V[1] = -1 ' Force Value of V[1] to be never equal to valid B
'ENDIF

' Now this is always true
'''IF B <> V[1] THEN

    ' VSOP basic values
    CLhvsop87a = LHVSOP87a(B)
    CRvsop87a = RVSOP87a(B)
    CBhvsop87a = BHVSOP87a(B)

'#DEBUG PRINT STR$(clhvsop87a)
'#DEBUG PRINT STR$(crvsop87a)
'#DEBUG PRINT STR$(cbhvsop87a)

    CDeltaPhi = DeltaPhiIAU(B)
    CDeltaEpsilon = DeltaEpsilonIAU(B)

    CLG = LH2LG(CLhvsop87a)

    CLtrue = LGVSOP2FK5(CLG) 'See Ltrue()

    CLambda = CLtrue + CDeltaPhi
    CLambda = CLambda - DMS2Unit(0, 0, 20.4898,%UnitDegree) / CRvsop87a ' See Lambda()

    CEpsilonMean = EpsilonMeanIAU(B)
    CEpsilon = CEpsilonMean + CDeltaEpsilon  ' See Epsilon()

    CBeta = BetaVSOP2FK5(BH2Beta(CBhvsop87a), B, CLG)

    CAlpha = LEB2Alpha(CLambda, CEpsilon, CBeta)

    CDelta = LEB2Delta(CLambda, CEpsilon, CBeta)

    CSunMeanLong = SunMeanLong(B)

    CEOT = Reduce180Signed(CSunMeanLong - 0.0057183 - CAlpha + CDeltaPhi * dCos(CEpsilon)) 'degrees
'    CEOT = CEOT - HMS2Unit(0,0,0.20,%UnitDegree) '0.2 sec adjustment to track MICA

    ' PUBLISHED
'''     V[1] = B ' Current Moment in time that is being calculated.
'''     V[2] = CEOT ' Equation of Time
'''     V[3] = CDelta ' Declination

    ' INTERNAL AND NOT AVAILABLE
'''     V[4] = CAlpha
'''     V[5] = CBeta

'''     V[6] = CEpsilon
'''     V[7] = CEpsilonMean
'''     V[8] = CDeltaPhi
'''     V[9] = CDeltaEpsilon
'''     V[10] = CLtrue

'''     V[11] = CLambda
'''     V[12] = CSunMeanLong
'''     V[13] = CLG
'''     V[14] = 0
'''     V[15] = CLhvsop87a

'''     V[16] = CRvsop87a
'''     V[17] = CBhvsop87a
'''ENDIF

'''IF SIndex < 1 || SIndex > 3 THEN SIndex = 1
'''FUNCTION = V[SIndex]
IF SIndex = 1 THEN
   FUNCTION = B
   EXIT FUNCTION
END IF
IF SIndex = 2 THEN
   FUNCTION = CEOT
   EXIT FUNCTION
END IF

FUNCTION = CDelta
END FUNCTION










'END #INCLUDE "astro.inc"        ' astronomical functions


'#INCLUDE "inttime.inc"

' #####################
' INTDeclMJDUnit
' #####################
'

FUNCTION INTDeclMJDUnit  (BYVAL MJD AS DOUBLE, BYVAL ToUnit AS LONG) AS DOUBLE
DIM  X AS EXT
X =Sun(Taum(JD2JDE(MJD2JD((MJD)))),3)
FUNCTION = DMS2ArcUnit(X,0.,0.,ToUnit)

END FUNCTION

' #######################
' INTDeclSDYUnit
' #######################
'

FUNCTION INTDeclSDYUnit(BYVAL SDY AS DOUBLE, BYVAL ToUnit AS LONG) AS DOUBLE
DIM RSDY   AS EXT
DIM X                AS EXT
DIM MJD                        AS EXT
DIM Feb01                                AS EXT
DIM Nov30                                          AS EXT
DIM MisJoin AS EXT

RSDY = Reduce1365((SDY))

MJD = YMD2MJD(%SDYIndexYear,1,(RSDY),0,0,0,0)

X = Sun(Taum(JD2JDE(MJD2JD((MJD)))),3)



MisJoin = SDYDeclMisJoin

Feb01 = YMD2SDY(2,1,0,0,0,0)

IF RSDY < Feb01 THEN ' January
    X = X + (Feb01-RSDY)/31 * MisJoin/2
END IF

Nov30 = YMD2SDY(11,30,0,0,0,0)
IF RSDY > Nov30 THEN  ' December
    X = X - (RSDY-Nov30)/32 * MisJoin /2
END IF


FUNCTION = DMS2ArcUnit(X,0.,0.,ToUnit)


END FUNCTION

' ####################
' INTEOTMJDUnit
' ####################
'

FUNCTION INTEOTMJDUnit(BYVAL MJD AS DOUBLE, BYVAL ToUnit AS LONG) AS DOUBLE

FUNCTION = DMS2Unit(Sun(Taum(JD2JDE(MJD2JD((MJD)))),2),0.,0.,ToUnit)
END FUNCTION

' #######################
' INTEOTSDYUnit
' #######################
'

FUNCTION INTEOTSDYUnit(BYVAL SDY AS DOUBLE, BYVAL ToUnit AS LONG) AS DOUBLE
DIM RSDY AS EXT
DIM X AS EXT
DIM MJD AS EXT
DIM Feb01         AS EXT
DIM Nov30                   AS EXT
DIM MisJoin AS EXT


RSDY = Reduce1365((SDY))

MJD = YMD2MJD(%SDYIndexYear,1,(RSDY),0,0,0,0)

' Blend values from previous year so that the cycle goes smoothly.
' this plan disturbs values through out the whole year.
X = Sun(Taum(JD2JDE(MJD2JD(MJD))),2) ' *(366.-RSDY)/365.+ Sun(Taum(JD2JDE(MJD2JD(MJD-365))),2)*(RSDY-1.)/365.

' Average values over 3 years so that the cycle goes smoothly.
' Still leaves a gap between end of Dec 31st and beginning of Jan 1
'X = (Sun(Taum(JD2JDE(MJD2JD(MJD))),2)+ Sun(Taum(JD2JDE(MJD2JD(MJD-365))),2)+ Sun(Taum(JD2JDE(MJD2JD(MJD+365))),2))/3.

' Apply modification to December and January's values to mate them up.

MisJoin = SDYEOTMisJoin

Feb01 = YMD2SDY(2,1,0,0,0,0)

IF RSDY < Feb01 THEN ' January
    X = X + (Feb01-RSDY)/31 * MisJoin/2
END IF

Nov30 = YMD2SDY(11,30,0,0,0,0)
IF RSDY > Nov30 THEN  ' December
    X = X - (RSDY-Nov30)/32 * MisJoin /2
END IF

FUNCTION = DMS2Unit(X,0.,0.,ToUnit)

END FUNCTION

'
' ########################
' INTUnit2Unit ()
' ########################
'
FUNCTION INTUnit2Unit(BYVAL Value AS DOUBLE, BYVAL FromUnit AS LONG, BYVAL ToUnit AS LONG)  AS DOUBLE

DIM  DValue  AS EXT



IF Value = %BadParameter THEN
   FUNCTION = %BadParameter
   EXIT FUNCTION
END IF

IF BadTimeArcUnit(ToUnit) THEN
   FUNCTION =  %BadParameter
   EXIT FUNCTION
END IF

IF BadTimeArcUnit(FromUnit) THEN
   FUNCTION = %BadParameter
   EXIT FUNCTION
END IF



' Convert to Degrees first
SELECT CASE FromUnit
    CASE %UnitDegree
        DValue = Value
    CASE %UnitArcMin
        DValue = Value/60.
    CASE %UnitArcSec
        DValue = Value/3600.
    CASE %UnitRadian
        DValue = Value*RADTODEG
    CASE %UnitDay
        DValue = Value*360.
    CASE %UnitHour
        DValue = Value*15.
    CASE %UnitTMin
        DValue = Value/4.
    CASE %UnitTSec
        DValue = Value/240.
END SELECT

' Convert to ToUnit
SELECT CASE ToUnit
    CASE %UnitDegree
        FUNCTION = DValue
        EXIT FUNCTION
    CASE %UnitArcMin
        FUNCTION = DValue * 60.
        EXIT FUNCTION
    CASE %UnitArcSec
        FUNCTION = DValue *3600.
        EXIT FUNCTION
    CASE %UnitRadian
        FUNCTION = DValue * DEGTORAD
        EXIT FUNCTION
    CASE %UnitDay
        FUNCTION = DValue /360.
        EXIT FUNCTION
    CASE %UnitHour
        FUNCTION = DValue/15.
        EXIT FUNCTION
    CASE %UnitTMin
        FUNCTION = DValue*4.
        EXIT FUNCTION
    CASE %UnitTSec
        FUNCTION = DValue*240.
        EXIT FUNCTION
END SELECT

FUNCTION = %BadParameter
END FUNCTION


' #######################
' NextX
' #######################
'
FUNCTION NextX(BYVAL SearchFor AS LONG, XMin AS DOUBLE, XMax AS DOUBLE, BYVAL YMin AS DOUBLE, BYVAL YGuess1 AS DOUBLE, BYVAL YGuess2 AS DOUBLE, BYVAL YMax AS DOUBLE, BYVAL E AS LONG) AS LONG

DIM XGuess1 AS EXT
DIM XGuess2 AS EXT

' Bad Parms check

IF SearchFor = %BadParameter THEN
   FUNCTION = %BadParameter
   EXIT FUNCTION
END IF
IF XMin = %BadParameter THEN
   FUNCTION = %BadParameter
   EXIT FUNCTION
END IF
IF XMax = %BadParameter THEN
   FUNCTION = %BadParameter
   EXIT FUNCTION
END IF
IF YMin = %BadParameter THEN
   FUNCTION = %BadParameter
   EXIT FUNCTION
END IF
IF YGuess1 = %BadParameter THEN
   FUNCTION = %BadParameter
   EXIT FUNCTION
END IF
IF YGuess2 = %BadParameter THEN
   FUNCTION = %BadParameter
   EXIT FUNCTION
END IF
IF YMax = %BadParameter THEN
   FUNCTION = %BadParameter
   EXIT FUNCTION
END IF
IF E = %BadParameter THEN
   FUNCTION = %BadParameter
   EXIT FUNCTION
END IF

'
' Parameter Checking
'
IF XMax <= XMin THEN
   FUNCTION = %BadParameter
   EXIT FUNCTION
END IF
IF XGuess1 <> 0. OR XGuess2 <> 0. THEN
    IF XMax <= XGuess2 THEN
       FUNCTION = %BadParameter
       EXIT FUNCTION
    END IF
    IF XGuess2 <= XGuess1 THEN
       FUNCTION = %BadParameter
        EXIT FUNCTION
    END IF
    IF XGuess1 <= XMin THEN
       FUNCTION = %BadParameter
        EXIT FUNCTION
    END IF
END IF
IF BadSearchFor(SearchFor) THEN
   FUNCTION = %BadParameter
   EXIT FUNCTION
END IF


'
' Check for found answer
'
IF SearchFor = %SearchZeroY  THEN
    IF ABS(YMax)<.9*10.^E THEN 'Found the best answer
'       #DEBUG PRINT YMax, YMin, YMax-YMin
        XMax = (XMax+XMin) / 2.
        XMin = XMax
        FUNCTION = %BinStatusFound
        EXIT FUNCTION
    END IF
ELSE
    IF (XMax-XMin)<.9*10.^E THEN 'Found the best answer
'       #DEBUG PRINT YMax, YMin, YMax-YMin
        XMax = (XMax+XMin) / 2.
        XMin = XMax
        FUNCTION = %BinStatusFound
        EXIT FUNCTION
    END IF
END IF

'
' Calculate New Guesses
'

'
' These two new X guesses MUST be generated identically by the calling program
XGuess1 = XMin+(XMax-XMin)*.5
XGuess2 = XMin+(XMax-XMin)*.51
SELECT CASE SearchFor
    CASE %SearchZeroY
        IF SGN(YMin) <> SGN(YGuess1) THEN
            XMax = XGuess1
        ELSE
            IF SGN(YGuess1) <> SGN(YGuess2) THEN
                XMin = XGuess1
                XMax = XGuess2
            ELSE
                IF SGN(YGuess2) <> SGN(YMax) THEN
                    XMin = XGuess2
                ELSE ' No solution possible
                    FUNCTION = %BinStatusNoSolution
                    EXIT FUNCTION
                END IF
            END IF
        END IF
    CASE %SearchMinY
        IF YGuess1 <= YGuess2 THEN
            XMax = XGuess2
        ELSE
            XMin = XGuess1
        END IF
    CASE %SearchMaxY
        IF YGuess1 <= YGuess2 THEN
            XMin = XGuess1
        ELSE
            XMax = XGuess2
        END IF
END SELECT
FUNCTION = %BinStatusSearching
END FUNCTION

'
' ########################
' INTFindAnnFeat ()
' ########################
'

FUNCTION INTFindAnnFeat(BYVAL Scale AS LONG, BYVAL Year AS LONG, BYVAL Feature AS LONG) AS DOUBLE
DIM  XMin AS DOUBLE
DIM XMax  AS DOUBLE
DIM XGuess1 AS DOUBLE
DIM XGuess2 AS DOUBLE
DIM YMin AS DOUBLE
DIM YMax  AS DOUBLE
DIM YGuess1 AS DOUBLE
DIM YGuess2 AS DOUBLE
DIM T  AS DOUBLE
DIM SearchFor AS LONG
DIM ToUnit AS LONG
DIM E AS LONG

DIM DA AS DOUBLE
DIM DB AS DOUBLE
DIM DDelta AS DOUBLE
DIM EOTA AS DOUBLE
DIM EOTB AS DOUBLE
DIM DeclA AS DOUBLE
DIM DeclB AS DOUBLE
DIM DTry AS DOUBLE
DIM EOTTry AS DOUBLE
DIM DeclTry AS DOUBLE
DIM DistAB AS DOUBLE
DIM DistATry AS DOUBLE
DIM DistBTry AS DOUBLE
DIM Swapped AS LONG
DIM Result AS LONG

ToUnit = %UnitDegree
E =-6

IF Feature = %FeatAnalemmaCrossing1 OR Feature = %FeatAnalemmaCrossing2 THEN
        IF Scale = %ScaleSDY THEN
            DA = YMD2SDY(10,1,0,0,0,0)
            DB = YMD2SDY(5,1,0,0,0,0)
        ELSE
            DA = YMD2MJD(Year,10,1,0,0,0,0)
            DB = YMD2MJD(Year,5,1,0,0,0,0)
        END IF
        IF Scale = %ScaleSDY THEN
            EOTA = INTEOTSDYUnit((DA),%UnitDegree)
            DeclA = INTDeclSDYUnit((DA),%UnitDegree)
            DeclB = INTDeclSDYUnit((DB),%UnitDegree)
            EOTB = INTEOTSDYUnit((DB),%UnitDegree)
        ELSE
            EOTA = INTEOTMJDUnit((DA),%UnitDegree)
            DeclA = INTDeclMJDUnit((DA),%UnitDegree)
            EOTB = INTEOTMJDUnit((DB),%UnitDegree)
            DeclB = INTDeclMJDUnit((DB),%UnitDegree)
        END IF
        DistAB = (EOTA-EOTB)^2+(DeclA-DeclB)^2
        DDelta = 15.
        DO
            Swapped = %FALSE
            DTry = DA+DDelta
            IF Scale = %ScaleSDY THEN
                EOTTry = INTEOTSDYUnit((DTry),%UnitDegree)
                DeclTry = INTDeclSDYUnit((DTry),%UnitDegree)
            ELSE
                EOTTry = INTEOTMJDUnit((DTry),%UnitDegree)
                DeclTry = INTDeclMJDUnit((DTry),%UnitDegree)
            END IF
            DistBTry =      (EOTTry-EOTB)^2+(DeclTry-DeclB)^2
            IF DistBTry < DistAB THEN
                DA = DTry
                EOTA = EOTTry
                DeclA = DeclTry
                DistAB = (EOTA-EOTB)^2+(DeclA-DeclB)^2
                Swapped = %TRUE
            ELSE
                DTry = DA-DDelta
                IF Scale = %ScaleSDY THEN
                    EOTTry = INTEOTSDYUnit((Dtry),%UnitDegree)
                    DeclTry = INTDeclSDYUnit((Dtry),%UnitDegree)
                ELSE
                    EOTTry = INTEOTMJDUnit((Dtry),%UnitDegree)
                    DeclTry = INTDeclMJDUnit((Dtry),%UnitDegree)
                END IF
                DistBTry =      (EOTTry-EOTB)^2+(DeclTry-DeclB)^2
                IF DistBTry < DistAB THEN
                    DA = DTry
                    EOTA = EOTTry
                    DeclA = DeclTry
                    DistAB = (EOTA-EOTB)^2+(DeclA-DeclB)^2
                    Swapped = %TRUE
                END IF
            END IF
            DTry = DB+DDelta
            IF Scale = %ScaleSDY THEN
                EOTTry = INTEOTSDYUnit((DTry),%UnitDegree)
                DeclTry = INTDeclSDYUnit((DTry),%UnitDegree)
            ELSE
                EOTTry = INTEOTMJDUnit((DTry),%UnitDegree)
                DeclTry = INTDeclMJDUnit((DTry),%UnitDegree)
            END IF
            DistATry =      (EOTTry-EOTA)^2+(DeclTry-DeclA)^2
            IF DistATry < DistAB THEN
                DB = DTry
                EOTB = EOTTry
                DeclB = DeclTry
                DistAB = (EOTA-EOTB)^2+(DeclA-DeclB)^2
                Swapped = %TRUE
            ELSE
                DTry = DB-DDelta
                IF Scale = %ScaleSDY THEN
                    EOTTry = INTEOTSDYUnit((DTry),%UnitDegree)
                    DeclTry = INTDeclSDYUnit((DTry),%UnitDegree)
                ELSE
                    EOTTry = INTEOTMJDUnit((DTry),%UnitDegree)
                    DeclTry = INTDeclMJDUnit((DTry),%UnitDegree)
                END IF
                DistATry =      (EOTTry-EOTA)^2+(DeclTry-DeclA)^2
                IF DistATry < DistAB THEN
                    DB = DTry
                    EOTB = EOTTry
                    DeclB = DeclTry
                    DistAB = (EOTA-EOTB)^2+(DeclA-DeclB)^2
                    Swapped = %TRUE
                END IF
            END IF
            IF ISFALSE Swapped THEN DDelta = DDelta/2.
        LOOP WHILE (DDelta)>.9*10.^E
        IF DA > DB THEN
            T = DA
            DA = DB
            DB = T
        END IF
        IF Feature = %FeatAnalemmaCrossing1 THEN
            FUNCTION = DA
            EXIT FUNCTION
        ELSE
            FUNCTION = DB
            EXIT FUNCTION
        END IF
'       #DEBUG PRINT DA, DB
END IF

DIM FuncToCalc AS LONG

SELECT  CASE Feature
    CASE %FeatEOTMax1
        ' Each Case does own setup of end points and type of search
        IF Scale = %ScaleSDY THEN
            XMin = YMD2SDY(5,1,0,0,0,0)
            XMax = YMD2SDY(5,31,0,0,0,0)
        ELSE
            XMin = YMD2MJD(Year,5,1,0,0,0,0)
            XMax = YMD2MJD(Year,5,31,0,0,0,0)
        END IF
        SearchFor = %SearchMaxY
        FuncToCalc = %FuncEOT
    CASE %FeatEOTMax2
        ' Each Case does own setup of end points and type of search
        IF Scale = %ScaleSDY THEN
            XMin = YMD2SDY(10,15,0,0,0,0)
            XMax = YMD2SDY(11,15,0,0,0,0)
        ELSE
            XMin = YMD2MJD(Year,10,15,0,0,0,0)
            XMax = YMD2MJD(Year,11,15,0,0,0,0)
        END IF
        SearchFor = %SearchMaxY
        FuncToCalc = %FuncEOT
    CASE %FeatEOTMin1
        ' Each Case does own setup of end points and type of search
        IF Scale = %ScaleSDY THEN
            XMin = YMD2SDY(1,15,0,0,0,0)
            XMax = YMD2SDY(2,15,0,0,0,0)
        ELSE
            XMin = YMD2MJD(Year,1,15,0,0,0,0)
            XMax = YMD2MJD(Year,2,15,0,0,0,0)
        END IF
        SearchFor = %SearchMinY
        FuncToCalc = %FuncEOT
    CASE %FeatEOTMin2
        ' Each Case does own setup of end points and type of search
        IF Scale = %ScaleSDY THEN
            XMin = YMD2SDY(7,15,0,0,0,0)
            XMax = YMD2SDY(8,15,0,0,0,0)
        ELSE
            XMin = YMD2MJD(Year,7,15,0,0,0,0)
            XMax = YMD2MJD(Year,8,15,0,0,0,0)
        END IF
        SearchFor = %SearchMinY
        FuncToCalc = %FuncEOT
    CASE %FeatEOTZero1
        ' Each Case does own setup of end points and type of search
        IF Scale = %ScaleSDY THEN
            XMin = YMD2SDY(4,1,0,0,0,0)
            XMax = YMD2SDY(4,30,0,0,0,0)
        ELSE
            XMin = YMD2MJD(Year,4,1,0,0,0,0)
            XMax = YMD2MJD(Year,4,30,0,0,0,0)
        END IF
        SearchFor = %SearchZeroY
        FuncToCalc = %FuncEOT
    CASE %FeatEOTZero2
        ' Each Case does own setup of end points and type of search
        IF Scale = %ScaleSDY THEN
            XMin = YMD2SDY(6,1,0,0,0,0)
            XMax = YMD2SDY(6,30,0,0,0,0)
        ELSE
            XMin = YMD2MJD(Year,6,1,0,0,0,0)
            XMax = YMD2MJD(Year,6,30,0,0,0,0)
        END IF
        SearchFor = %SearchZeroY
        FuncToCalc = %FuncEOT
    CASE %FeatEOTZero3
        ' Each Case does own setup of end points and type of search
        IF Scale = %ScaleSDY THEN
            XMin = YMD2SDY(8,15,0,0,0,0)
            XMax = YMD2SDY(9,15,0,0,0,0)
        ELSE
            XMin = YMD2MJD(Year,8,15,0,0,0,0)
            XMax = YMD2MJD(Year,9,15,0,0,0,0)
        END IF
        SearchFor = %SearchZeroY
        FuncToCalc = %FuncEOT
    CASE %FeatEOTZero4
        ' Each Case does own setup of end points and type of search
        IF Scale = %ScaleSDY THEN
            XMin = YMD2SDY(12,1,0,0,0,0)
            XMax = YMD2SDY(1,30,0,0,0,0)+365
        ELSE
            XMin = YMD2MJD(Year,12,1,0,0,0,0)
            XMax = YMD2MJD(Year+1,1,30,0,0,0,0)
        END IF
        SearchFor = %SearchZeroY
        FuncToCalc = %FuncEOT
    CASE %FeatDeclMax
        ' Each Case does own setup of end points and type of search
        IF Scale = %ScaleSDY THEN
            XMin = YMD2SDY(6,1,0,0,0,0)
            XMax = YMD2SDY(7,7,0,0,0,0)
        ELSE
            XMin = YMD2MJD(Year,6,1,0,0,0,0)
            XMax = YMD2MJD(Year,7,7,0,0,0,0)
        END IF
        SearchFor = %SearchMaxY
        FuncToCalc = %FuncDecl
    CASE %FeatDeclMin
        ' Each Case does own setup of end points and type of search
        IF Scale = %ScaleSDY THEN
            XMin = YMD2SDY(12,1,0,0,0,0)
            XMax = YMD2SDY(1,15,0,0,0,0)+365
        ELSE
            XMin = YMD2MJD(Year,12,1,0,0,0,0)
            XMax = YMD2MJD(Year+1,1,15,0,0,0,0)
        END IF
        SearchFor = %SearchMinY
        FuncToCalc = %FuncDecl
    CASE %FeatDeclZero1
        ' Each Case does own setup of end points and type of search
        IF Scale = %ScaleSDY THEN
            XMin = YMD2SDY(1,1,0,0,0,0)
            XMax = YMD2SDY(6,1,0,0,0,0)
        ELSE
            XMin = YMD2MJD(Year,1,1,0,0,0,0)
            XMax = YMD2MJD(Year,6,1,0,0,0,0)
        END IF
        SearchFor = %SearchZeroY
        FuncToCalc = %FuncDecl
    CASE %FeatDeclZero2
        ' Each Case does own setup of end points and type of search
        IF Scale = %ScaleSDY THEN
            XMin = YMD2SDY(6,1,0,0,0,0)
            XMax = YMD2SDY(12,1,0,0,0,0)
        ELSE
            XMin = YMD2MJD(Year,6,1,0,0,0,0)
            XMax = YMD2MJD(Year,12,1,0,0,0,0)
        END IF
        SearchFor = %SearchZeroY
        FuncToCalc = %FuncDecl

    CASE %FeatEquinox1
        ' Each Case does own setup of end points and type of search
        IF Scale = %ScaleSDY THEN
            XMin = YMD2SDY(3,15,0,0,0,0)
            XMax = YMD2SDY(3,30,0,0,0,0)
        ELSE
            XMin = YMD2MJD(Year,3,15,0,0,0,0)
            XMax = YMD2MJD(Year,3,30,0,0,0,0)
        END IF
        SearchFor = %SearchZeroY
        FuncToCalc = %FuncEquinox1

    CASE %FeatSolstice1
        ' Each Case does own setup of end points and type of search
        IF Scale = %ScaleSDY THEN
            XMin = YMD2SDY(6,15,0,0,0,0)
            XMax = YMD2SDY(6,30,0,0,0,0)
        ELSE
            XMin = YMD2MJD(Year,6,15,0,0,0,0)
            XMax = YMD2MJD(Year,6,30,0,0,0,0)
        END IF
        SearchFor = %SearchZeroY
        FuncToCalc = %FuncSolstice1

    CASE %FeatEquinox2
        ' Each Case does own setup of end points and type of search
        IF Scale = %ScaleSDY THEN
            XMin = YMD2SDY(9,15,0,0,0,0)
            XMax = YMD2SDY(9,30,0,0,0,0)
        ELSE
            XMin = YMD2MJD(Year,9,15,0,0,0,0)
            XMax = YMD2MJD(Year,9,30,0,0,0,0)
        END IF
        SearchFor = %SearchZeroY
        FuncToCalc = %FuncEquinox2

    CASE %FeatSolstice2
        ' Each Case does own setup of end points and type of search
        IF Scale = %ScaleSDY THEN
            XMin = YMD2SDY(12,15,0,0,0,0)
            XMax = YMD2SDY(12,30,0,0,0,0)
        ELSE
            XMin = YMD2MJD(Year,12,15,0,0,0,0)
            XMax = YMD2MJD(Year,12,30,0,0,0,0)
        END IF
        SearchFor = %SearchZeroY
        FuncToCalc = %FuncSolstice2

END SELECT
' Check units for this FuncToCalc

IF FuncToCalc = %FuncEOT THEN
    IF BadTimeArcUnit(ToUnit) THEN
       FUNCTION = %BadParameter
       EXIT FUNCTION
    END IF
ELSE
    IF BadArcUnit(ToUnit) THEN
       FUNCTION = %BadParameter
       EXIT FUNCTION
    END IF

END IF

DIM SolsticeDiff AS LONG

DO
' These two new X guesses MUST be generated identically here and by the NextX() Function
    XGuess1 = XMin+(XMax-XMin)*.5
    XGuess2 = XMin+(XMax-XMin)*.51
    SELECT CASE FuncToCalc
        CASE %FuncEOT
           IF Scale = %ScaleSDY THEN
                YMin = INTEOTSDYUnit((XMin),ToUnit)
                YGuess1 = INTEOTSDYUnit((XGuess1),ToUnit)
                YGuess2 = INTEOTSDYUnit((XGuess2),ToUnit)
                YMax = INTEOTSDYUnit((XMax),ToUnit)
           ELSE
                YMin = INTEOTMJDUnit((XMin),ToUnit)
                YGuess1 = INTEOTMJDUnit((XGuess1),ToUnit)
                YGuess2 = INTEOTMJDUnit((XGuess2),ToUnit)
                YMax = INTEOTMJDUnit((XMax),ToUnit)
           END IF
        CASE %FuncDecl
             IF  Scale = %ScaleSDY THEN
                YMin = INTDeclSDYUnit((XMin),ToUnit)
                YGuess1 = INTDeclSDYUnit((XGuess1),ToUnit)
                YGuess2 = INTDeclSDYUnit((XGuess2),ToUnit)
                YMax = INTDeclSDYUnit((XMax),ToUnit)
             ELSE
                YMin = INTDeclMJDUnit((XMin),ToUnit)
                YGuess1 = INTDeclMJDUnit((XGuess1),ToUnit)
                YGuess2 = INTDeclMJDUnit((XGuess2),ToUnit)
                YMax = INTDeclMJDUnit((XMax),ToUnit)
             END IF
        CASE %FuncSolstice1, %FuncSolstice2, %FuncEquinox1, %FuncEquinox2
             SELECT CASE FuncToCalc
                    CASE %FuncEquinox1
                         SolsticeDiff = 0
                    CASE %FuncSolstice1
                         SolsticeDiff = 90
                    CASE %FuncEquinox2
                         SolsticeDiff = 180
                    CASE %FuncSolstice2
                         SolsticeDiff = 270
             END SELECT
             IF  Scale = %ScaleSDY THEN
                YMin = reduce180signed(lambda(taum(jd2jde(mjd2jd(SDY2MJD(%SDYindexyear,XMin))))) - SolsticeDiff)
                YGuess1 = reduce180signed(lambda(taum(jd2jde(mjd2jd(SDY2MJD(%SDYindexyear,XGuess1))))) - SolsticeDiff)
                YGuess2 = reduce180signed(lambda(taum(jd2jde(mjd2jd(SDY2MJD(%SDYindexyear,XGuess2))))) - SolsticeDiff)
                YMax = reduce180signed(lambda(taum(jd2jde(mjd2jd(SDY2MJD(%SDYindexyear,XMax))))) - SolsticeDiff)
             ELSE
                YMin = reduce180signed(lambda(taum(jd2jde(mjd2jd(XMin)))) - SolsticeDiff)
                YGuess1 = reduce180signed(lambda(taum(jd2jde(mjd2jd(XGuess1)))) - SolsticeDiff)
                YGuess2 = reduce180signed(lambda(taum(jd2jde(mjd2jd(XGuess2)))) - SolsticeDiff)
                YMax = reduce180signed(lambda(taum(jd2jde(mjd2jd(XMax)))) - SolsticeDiff)
' These lines for testing give result in TD cf. p182 Meeus +/-10sec VSOP87
'                YMin = reduce180signed(lambda(taum(mjd2jd(XMin))) - SolsticeDiff)
'                YGuess1 = reduce180signed(lambda(taum(mjd2jd(XGuess1))) - SolsticeDiff)
'                YGuess2 = reduce180signed(lambda(taum(mjd2jd(XGuess2))) - SolsticeDiff)
'                YMax = reduce180signed(lambda(taum(mjd2jd(XMax))) - SolsticeDiff)
             END IF
    END SELECT
    ' NOTE WELL: Two parameters MUST called by reference otherwise infinite loop.
    Result = NextX(SearchFor, XMin, XMax, YMin, YGuess1, YGuess2, YMax, E)
    IF XMin = %BadParameter THEN
       FUNCTION = %BadParameter ' Check if NextX was fed a wrong parameter
       EXIT FUNCTION
    END IF
    IF XMax = %BadParameter THEN
       FUNCTION = %BadParameter
       EXIT FUNCTION
    END IF
LOOP WHILE Result = %BinStatusSearching

IF Result = %BinStatusFound THEN
    FUNCTION = XMin
    EXIT FUNCTION
ELSE
    FUNCTION = %BinStatusNoSolution
    EXIT FUNCTION
END IF

END FUNCTION

FUNCTION INTH(BYVAL SunType AS LONG, BYVAL Time AS DOUBLE, BYVAL Longitude AS DOUBLE) AS DOUBLE

DIM H AS DOUBLE

H = ROUND(INTUnit2Unit(Time-INT(Time),%UnitDay,%UnitDegree),8)
SELECT CASE SunType
    CASE %SunMeanLocal
        FUNCTION = Reduce180Signed(H-180.) + V(8)*(0.75*dSin(h+56)+0.25*dSin(h*3.14159))
        EXIT FUNCTION
    CASE %SunMeanAtLong
        FUNCTION = Reduce180Signed(H-180.- Longitude)+ V(8)*(0.75*dSin(h+23)+0.25*dSin(h*5.7897980))
        EXIT FUNCTION
    CASE %SunSDYAppAtLong
        FUNCTION = Reduce180Signed(H-180.- Longitude + INTEOTSDYUnit(Time,%UnitDegree))
        EXIT FUNCTION
    CASE %SunMJDAppAtLong
        FUNCTION = Reduce180Signed(H-180.- Longitude + INTEOTMJDUnit(Time,%UnitDegree))
        EXIT FUNCTION
END SELECT

END FUNCTION

'
' #######################
' INTHPD2Azi
' #######################
'
FUNCTION INTHPD2Azi(BYVAL H AS DOUBLE, BYVAL Phi AS DOUBLE, BYVAL  Decl AS DOUBLE) AS DOUBLE
DIM A AS EXT
DIM B AS EXT
DIM Azi AS EXT
DIM xH AS EXT

IF H = 0 AND Phi = Decl THEN  ' Sun is straight overhead and Azi is undefined
   IF Phi > 0 THEN
        FUNCTION = 180
   ELSE
        FUNCTION = 0
   END IF
   EXIT FUNCTION
END IF

xH = Reduce180Signed((H))

A = dSin(xH)
B = dCos(xH)*dSin((Phi))- dTan((Decl))*dCos(Phi)

Azi = dAtan(A/B)+ 180.
IF B < 0. THEN
    Azi = Azi + 180.
END IF

DIM result AS DOUBLE
result = Reduce360(Azi)
FUNCTION = result

END FUNCTION

' ####################
' INTAziMJDUnit
' ####################
'


FUNCTION INTAziMJDUnit ALIAS "INTAziMJDUnit" (BYVAL MJD AS DOUBLE, BYVAL Longitude AS DOUBLE, BYVAL Phi AS DOUBLE, BYVAL ToUnit AS LONG) AS DOUBLE
FUNCTION = INTUnit2Unit(INTHPD2Azi(INTH(%SunMJDAppAtLong,MJD,Longitude), Phi, INTDeclMJDUnit(MJD,%UnitDegree)),%UnitDegree,ToUnit)
END FUNCTION
' ####################
' sdAziSDYUnit
' ####################
'

FUNCTION INTAziSDYUnit(BYVAL SDY AS DOUBLE, BYVAL Longitude AS DOUBLE, BYVAL Phi AS DOUBLE, BYVAL ToUnit AS LONG) AS DOUBLE
DIM xSDY AS EXT

xSDY = Reduce1365((SDY))

FUNCTION = INTUnit2Unit(INTHPD2Azi(INTH(%SunSDYAppAtLong,(xSDY),(Longitude)), (Phi), INTDeclSDYUnit((xSDY),%UnitDegree)),%Unitdegree,ToUnit)

END FUNCTION

'
' #######################
' INTHPD2Alt
' #######################
'
FUNCTION INTHPD2Alt(BYVAL H AS DOUBLE, BYVAL Phi AS DOUBLE, BYVAL  Decl AS DOUBLE) AS DOUBLE
DIM A AS EXT
DIM B AS EXT
DIM Alt AS EXT
DIM xH AS EXT

xH = Reduce180Signed((H))

Alt = dASin(dSin(Phi)*dSin(Decl)+dCos(Phi)*dCos(Decl)*dCos(xH))

'Function = Reduce360(Alt) ' This was a misconception removed!
FUNCTION = Alt

END FUNCTION

' ####################
' INTAltMJDUnit
' ####################
'

FUNCTION INTAltMJDUnit(BYVAL MJD AS DOUBLE, BYVAL Longitude AS DOUBLE, BYVAL Phi AS DOUBLE, BYVAL ToUnit AS LONG) AS DOUBLE
FUNCTION = INTUnit2Unit(INTHPD2Alt(INTH(%SunMJDAppAtLong,MJD,Longitude), Phi, INTDeclMJDUnit(MJD,%UnitDegree)),%UnitDegree,ToUnit)
END FUNCTION
' ####################
' sdAziSDYUnit
' ####################
'

FUNCTION INTAltSDYUnit(BYVAL SDY AS DOUBLE, BYVAL Longitude AS DOUBLE, BYVAL Phi AS DOUBLE, BYVAL ToUnit AS LONG) AS DOUBLE
DIM xSDY AS EXT

xSDY = Reduce1365((SDY))

FUNCTION = INTUnit2Unit(INTHPD2Alt(INTH(%SunSDYAppAtLong,(xSDY),(Longitude)), (Phi), INTDeclSDYUnit((xSDY),%UnitDegree)),%Unitdegree,ToUnit)

END FUNCTION


'END #INCLUDE "inttime.inc"

FUNCTION LoadKeyFile() AS LONG

DIM LastSerial AS STRING
DIM LastDateFlash AS STRING
DIM KeyFile AS STRING
DIM ret AS LONG
DIM V(%RegNoofVs)

LastDateFlash = GetRegValue(%HKEY_LOCAL_MACHINE, "SOFTWARE\"+$RegCompany+"\"+$RegProduct+"\"+$RegVersion, $RegValueLastDateFlash)
LastDateFlash  = RC4Like(HexS2S(LastDateFlash),StaticKey())


LastSerial = GetRegValue(%HKEY_LOCAL_MACHINE, "SOFTWARE\"+$RegCompany+"\"+$RegProduct+"\"+$RegVersion, $RegValueSerial)
LastSerial  = RC4Like(HexS2S(LastSerial),StaticKey())

keyfile = GetRegValue(%HKEY_LOCAL_MACHINE, "SOFTWARE\"+$RegCompany+"\"+$RegProduct+"\"+$RegVersion, $RegValueKeyFile)

keyfile = xgets($KeyfileTagKeyFile,Keyfile)
KeyFile = RC4Like(HexS2S(keyfile),StaticKey())

'MsgBox keyfile

Serial  = xgets($KeyFileTagSerial,KeyFile)
IF Serial = "" THEN 'most likely corrupted KeyFile
    MSGBOX "Load Failure 181"+$cr+$lf+$cr+$lf _
        , %MB_OK, $DialogTitle
    FUNCTION = %FALSE
    EXIT FUNCTION
END IF

KeyFileLevel = xgetv($KeyFileTagKeyFileLevel,KeyFile)
IF KeyFileLevel = 0 THEN 'most likely corrupted KeyFile
    MSGBOX "Load Failure 182"+$cr+$lf+$cr+$lf _
        , %MB_OK,  $DialogTitle
    FUNCTION = %FALSE
    EXIT FUNCTION
END IF

IF KeyFileLevel <> 1 THEN 'most likely keyfile from future version
    MSGBOX "Load Failure 200"+$cr+$lf+$cr+$lf _
        , %MB_OK,  $DialogTitle
    FUNCTION = %FALSE
    EXIT FUNCTION
END IF


KeyType  = xgets($KeyFileTagKeyType,KeyFile)
IF KeyType = "" THEN 'most likely corrupted KeyFile
    MSGBOX "Load Failure 182/1"+$cr+$lf+$cr+$lf _
        , %MB_OK, $DialogTitle
    FUNCTION = %FALSE
    EXIT FUNCTION
END IF

Company  = xgets($KeyFileTagCompany,KeyFile)

Email  = xgets($KeyFileTagEmail,KeyFile)
IF Email = "" THEN 'most likely corrupted KeyFile
    MSGBOX "Load Failure 182/3"+$cr+$lf+$cr+$lf _
        , %MB_OK, $DialogTitle
    FUNCTION = %FALSE
    EXIT FUNCTION
END IF

TimeStamp  = xgets($KeyFileTagTimeStamp,KeyFile)
IF TimeStamp = "" THEN 'most likely corrupted KeyFile
    MSGBOX "Load Failure 183"+$cr+$lf+$cr+$lf _
        , %MB_OK, $DialogTitle
    FUNCTION = %FALSE
    EXIT FUNCTION
END IF

User  = xgets($KeyFileTagUser,KeyFile)
IF User = "" THEN 'most likely corrupted KeyFile
    MSGBOX "Load Failure 184"+$cr+$lf+$cr+$lf _
        , %MB_OK, $DialogTitle
    FUNCTION = %FALSE
    EXIT FUNCTION
END IF

DIM i AS LONG


FOR i = 1 TO %RegNoofVs
    V(i) = xgetv($KeyFileTagV+FORMAT$(i),KeyFile)
NEXT

IF LastSerial <> Serial AND LastSerial<>"" THEN
   MSGBOX "Thankyou - you have just installed a new version or license of The Sun API." _
        , %MB_OK, $DialogTitle
   ret = SetRegValue(%HKEY_LOCAL_MACHINE, "SOFTWARE\"+$RegCompany+"\"+$RegProduct+"\"+$RegVersion, $RegValueSerial, S2HexS(RC4Like(Serial,StaticKey())))
   LastDateFlash = ""
END IF

IF LastDateFlash <> DATE$ THEN
    ret = sdinfo()
    ret = SetRegValue(%HKEY_LOCAL_MACHINE, "SOFTWARE\"+$RegCompany+"\"+$RegProduct+"\"+$RegVersion, $RegValueLastDateFlash, S2HexS(RC4Like(DATE$,StaticKey())))
END IF

FUNCTION = %TRUE



END FUNCTION



SUB InitGlobals
EarliestMJD = VAL($EarliestMJD)
LatestMJD = VAL($LatestMJD)
SDYEOTMisJoin = 0
SDYEOTMisJoin =INTEOTSDYUnit(365.999999999,%UnitDegree)-INTEOTSDYUnit(1,%UnitDegree)
SDYDeclMisJoin = 0
SDYDeclMisJoin =INTDeclSDYUnit(365.999999999,%UnitDegree)-INTDeclSDYUnit(1,%UnitDegree)




END SUB



'*************************************************************
'*************************************************************
'*************************************************************
'START OF EXPORTS
'*************************************************************
'*************************************************************
'*************************************************************

' WARNING all STRING type exported parameters must be declared
'   in PB as "... as Asciiz" and be
'   matched In VB Function Declare with "ByVal ... As String"


' #######################
' sdVersion
' #######################
'

FUNCTION sdVersion ALIAS "sdVersion" () EXPORT AS STRING
        FUNCTION = $version+" "+keytype
END FUNCTION


' #######################
' sdInfo
' #######################
'

FUNCTION sdInfo ALIAS "sdInfo" () EXPORT AS LONG

DIM typemsg AS STRING

IF keytype = $KeyTypeFreeware THEN
    typemsg = _
            "Freeware Accuracy:"+$cr+$lf+ _
            "   Equation of Time +/- 1 minute"+$cr+$lf+ _
            "   Declination      +/- 1 degree"+$cr+$lf+ _
            "   Hour Angle       +/- 1 degree"+$cr+$lf+$cr+$lf+ _
            "Professional Accuracy Version available at"+$cr+$lf+ _
            "   www.sunlit-design.com"
ELSE
   typemsg = _
            "Serial No: "+serial+$cr+$lf+$cr+$lf+ _
            "LICENSED TO: "+$cr+$lf+ _
            "   User: "+user+$cr+$lf+ _
            "   Company: "+Company+$cr+$lf+ _
            "   Email: "+email

END IF

        MSGBOX sdversion()+$cr+$lf+$cr+$lf+ _
            "*** "+keytype+" Version ***"+$cr+$lf+$cr+$lf+ _
            typemsg +$cr+$lf+$cr+$lf+ _
            $DialogCopyRight,_
            %MB_OK, $DialogTitle
END FUNCTION

' #######################
' sdMJDx
' #######################
'

FUNCTION sdMJDx ALIAS "sdMJDx" (BYVAL Year AS LONG, BYVAL Month AS LONG, BYVAL  Day AS DOUBLE, BYVAL Hour AS DOUBLE, BYVAL Minute AS DOUBLE, BYVAL Sec AS DOUBLE, BYVAL TimeZone AS DOUBLE) EXPORT AS STRING

DIM X AS EXT

IF Year = %BadParameter THEN
   FUNCTION = xerror(%BadParameter)
   EXIT FUNCTION
END IF
IF Month = %BadParameter THEN
   FUNCTION = xerror(%BadParameter)
   EXIT FUNCTION
END IF
IF Day = %BadParameter THEN
   FUNCTION = xerror(%BadParameter)
   EXIT FUNCTION
END IF
IF Hour = %BadParameter THEN
   FUNCTION = xerror(%BadParameter)
   EXIT FUNCTION
END IF
IF Minute = %BadParameter THEN
   FUNCTION = xerror(%BadParameter)
   EXIT FUNCTION
END IF
IF Sec = %BadParameter THEN
   FUNCTION = xerror(%BadParameter)
   EXIT FUNCTION
END IF
IF TimeZone = %BadParameter THEN
   FUNCTION = xerror(%BadParameter)
   EXIT FUNCTION
END IF

IF BadDate(Year,Month,1) THEN
   FUNCTION = xerror(%BadParameter)
   EXIT FUNCTION
END IF

IF BadTimeZone(TimeZone) THEN
   FUNCTION = xerror(%BadParameter)
   EXIT FUNCTION
END IF

IF Year=1582 AND Month=10 AND 5 <= Day AND Day <= 14 THEN 'Gregorian Reform
   FUNCTION = xerror(%BadParameter)
   EXIT FUNCTION
END IF

DIM GregDay AS DOUBLE

GregDay = Day

IF Year=1582 AND Month=10 AND GregDay > 14 THEN 'Gregorian Reform
   GregDay = GregDay - 10
END IF

X = YMD2MJD(Year, Month , GregDay , Hour, Minute, Sec, TimeZone)

IF BadMJD(X) THEN
   FUNCTION = xerror(%BadParameter)
   EXIT FUNCTION
END IF

FUNCTION = STR$(X)+xwrap($TagMJD,STR$(X))+xwrap($TagTimeZone,STR$(TimeZone))

END FUNCTION


' #####################
' sdSDYx
' #####################
'

FUNCTION sdSDYx ALIAS "sdSDYx" (BYVAL Month AS LONG, BYVAL Day AS DOUBLE, BYVAL Hour AS DOUBLE, BYVAL Minute AS DOUBLE, BYVAL Sec AS DOUBLE, BYVAL TimeZone AS DOUBLE) EXPORT AS STRING

IF Month = %BadParameter THEN
   FUNCTION = xerror(%BadParameter)
   EXIT FUNCTION
END IF
IF Day = %BadParameter THEN
   FUNCTION = xerror(%BadParameter)
   EXIT FUNCTION
END IF
IF Hour = %BadParameter THEN
   FUNCTION = xerror(%BadParameter)
   EXIT FUNCTION
END IF
IF Minute = %BadParameter THEN
   FUNCTION = xerror(%BadParameter)
   EXIT FUNCTION
END IF
IF Sec = %BadParameter THEN
   FUNCTION = xerror(%BadParameter)
   EXIT FUNCTION
END IF
IF TimeZone = %BadParameter THEN
   FUNCTION = xerror(%BadParameter)
   EXIT FUNCTION
END IF

IF BadDate(%SDYindexyear,Month,1) THEN
   FUNCTION = xerror(%BadParameter)
   EXIT FUNCTION
END IF

IF BadTimeZone(TimeZone) THEN
   FUNCTION = xerror(%BadParameter)
   EXIT FUNCTION
END IF

DIM value AS STRING

value = STR$(YMD2SDY(Month , Day , Hour, Minute, Sec, TimeZone))

FUNCTION = value+" "+xwrap($TagSDY,value)+xwrap($TagTimeZone,STR$(TimeZone))

END FUNCTION

'
' ########################
' sdCTDx ()
' ########################
'
FUNCTION sdCTDx ALIAS "sdCTDx" (BYVAL Hour AS DOUBLE, BYVAL Minute AS DOUBLE,BYVAL Sec AS DOUBLE) EXPORT AS STRING

IF Hour = %BadParameter THEN
   FUNCTION = xerror(%BadParameter)
   EXIT FUNCTION
END IF
IF Minute = %BadParameter THEN
   FUNCTION = xerror(%BadParameter)
   EXIT FUNCTION
END IF
IF Sec = %BadParameter THEN
   FUNCTION = xerror(%BadParameter)
   EXIT FUNCTION
END IF

DIM value AS STRING

value = STR$(HMS2CTD((Hour),(Minute),(Sec)))
FUNCTION = value+xwrap($TagCTD,value)

END FUNCTION

'
' ########################
' sdUTDx ()
' ########################
'
FUNCTION sdUTDx ALIAS "sdUTDx" (BYVAL Hour AS DOUBLE, BYVAL Minute AS DOUBLE,BYVAL Sec AS DOUBLE, BYVAL TimeZone AS DOUBLE) EXPORT AS STRING

IF Hour = %BadParameter THEN
   FUNCTION = xerror(%BadParameter)
   EXIT FUNCTION
END IF
IF Minute = %BadParameter THEN
   FUNCTION = xerror(%BadParameter)
   EXIT FUNCTION
END IF
IF Sec = %BadParameter THEN
   FUNCTION = xerror(%BadParameter)
   EXIT FUNCTION
END IF

IF TimeZone = %BadParameter THEN
   FUNCTION = xerror(%BadParameter)
   EXIT FUNCTION
END IF

DIM value AS STRING

value = STR$(HMST2UTD((Hour),(Minute),(Sec),(TimeZone)) )
FUNCTION = value+xwrap($TagUTD,value)+xwrap($TagTimeZone,STR$(TimeZone))

END FUNCTION


' #####################
' sdxDecl
' #####################
'

FUNCTION sdxDecl ALIAS "sdxDecl" (Timex AS ASCIIZ) EXPORT AS DOUBLE

IF xgetv($TagError,Timex) = %BadParameter THEN
    FUNCTION = %BadParameter
    EXIT FUNCTION
END IF

DIM xid AS STRING
DIM xvalue AS EXT
xvalue = xUnWrapTime(xid,Timex)
IF xvalue = %BadParameter THEN
   FUNCTION = %BadParameter
   EXIT FUNCTION
END IF


SELECT CASE xid
   CASE $TagMJD
        IF BadMJD(xvalue) THEN
           FUNCTION = %BadParameter
           EXIT FUNCTION
        END IF
        FUNCTION = INTDeclMJDUnit(xvalue,%UnitDegree)
   CASE $TagSDY
        IF BadSDY(xvalue) THEN
           FUNCTION = %BadParameter
           EXIT FUNCTION
        END IF
        FUNCTION = INTDeclSDYUnit(xvalue,%UnitDegree)
    CASE ELSE
        FUNCTION = %BadParameter
    EXIT FUNCTION
END SELECT

END FUNCTION


' #####################
' sdxEOT
' #####################
'

FUNCTION sdxEOT ALIAS "sdxEOT" (timex AS ASCIIZ) EXPORT AS DOUBLE
DIM  X AS EXT

IF xgetv($TagError,Timex) THEN
    FUNCTION = %BadParameter
    EXIT FUNCTION
END IF

DIM xid AS STRING
DIM xvalue AS EXT
xvalue = xUnWrapTime(xid,Timex)
IF xvalue = %BadParameter THEN
   FUNCTION = %BadParameter
   EXIT FUNCTION
END IF


SELECT CASE xid
   CASE $TagMJD
        IF BadMJD(xvalue) THEN
           FUNCTION = %BadParameter
           EXIT FUNCTION
        END IF
        FUNCTION = INTEOTMJDUnit(xvalue,%UnitDegree)
   CASE $TagSDY
        IF BadSDY(xvalue) THEN
           FUNCTION = %BadParameter
           EXIT FUNCTION
        END IF
        FUNCTION = INTEOTSDYUnit(xvalue,%UnitDegree)
    CASE ELSE
        FUNCTION = %BadParameter
    EXIT FUNCTION
END SELECT

END FUNCTION

'
' ########################
' sdxH ()
' ########################
'

FUNCTION sdxH ALIAS "sdxH" (timex AS ASCIIZ,BYVAL Longitude AS DOUBLE) EXPORT AS DOUBLE

DIM H AS DOUBLE

IF xgetv($TagError,Timex) THEN
   FUNCTION = %BadParameter
   EXIT FUNCTION
END IF
IF Longitude = %BadParameter THEN
   FUNCTION = %BadParameter
   EXIT FUNCTION
END IF

DIM xid AS STRING
DIM xvalue AS EXT
xvalue = xUnWrapTime(xid,Timex)
IF xvalue = %BadParameter THEN
   FUNCTION = %BadParameter
   EXIT FUNCTION
END IF


IF xid = $TagMJD OR xid = $TagSDY OR xid = $TagUTD THEN
    IF BadLongitude(Longitude) THEN
       FUNCTION = %BadParameter
       EXIT FUNCTION
    END IF
END IF

SELECT CASE xid
   CASE $TagCTD
        IF BadCTD(xvalue) THEN
           FUNCTION = %BadParameter
           EXIT FUNCTION
        END IF
        FUNCTION = INTH(%SunMeanLocal,xvalue,Longitude)
        EXIT FUNCTION
   CASE $TagUTD
        IF BadUTD(xvalue) THEN
           FUNCTION = %BadParameter
           EXIT FUNCTION
        END IF
        FUNCTION = INTH(%SunMeanAtLong,xvalue,Longitude)
        EXIT FUNCTION
   CASE $TagSDY
        IF BadSDY(xvalue) THEN
           FUNCTION = %BadParameter
           EXIT FUNCTION
        END IF
        FUNCTION = INTH(%SunSDYAppAtLong,xvalue,Longitude)
        EXIT FUNCTION
   CASE $TagMJD
        IF BadMJD(xvalue) THEN
           FUNCTION = %BadParameter
           EXIT FUNCTION
        END IF
        FUNCTION = INTH(%SunMJDAppAtLong,xvalue,Longitude)
        EXIT FUNCTION
   CASE ELSE
        FUNCTION = %BadParameter
        EXIT FUNCTION
END SELECT

END FUNCTION

' ####################
' sdxAzi
' ####################
'

FUNCTION sdxAzi ALIAS "sdxAzi" (timex AS ASCIIZ,BYVAL Longitude AS DOUBLE, BYVAL Latitude AS DOUBLE) EXPORT AS DOUBLE

IF xgetv($TagError,Timex) THEN
   FUNCTION = %BadParameter
   EXIT FUNCTION
END IF
IF Longitude = %BadParameter THEN
   FUNCTION = %BadParameter
   EXIT FUNCTION
END IF
IF Latitude = %BadParameter THEN
   FUNCTION = %BadParameter
   EXIT FUNCTION
END IF

DIM xid AS STRING
DIM xvalue AS EXT
xvalue = xUnWrapTime(xid,Timex)
IF xvalue = %BadParameter THEN
   FUNCTION = %BadParameter
   EXIT FUNCTION
END IF

IF BadLatitude(Latitude) THEN
   FUNCTION = %BadParameter
   EXIT FUNCTION
END IF
IF BadLongitude(Longitude) THEN
   FUNCTION = %BadParameter
   EXIT FUNCTION
END IF

SELECT CASE xid
       CASE $TagMJD
            IF BadMJD(xvalue) THEN
               FUNCTION = %BadParameter
               EXIT FUNCTION
            END IF
            FUNCTION = INTAziMJDUnit(xvalue, Longitude , Latitude, %Unitdegree)
            EXIT FUNCTION
       CASE $TagSDY
            IF BadSDY(xvalue) THEN
               FUNCTION = %BadParameter
               EXIT FUNCTION
            END IF
            FUNCTION = INTAziSDYUnit(xvalue, Longitude , Latitude, %Unitdegree)
            EXIT FUNCTION
       CASE ELSE
            FUNCTION = %BadParameter
            EXIT FUNCTION
END SELECT
END FUNCTION

' ####################
' sdxAlt
' ####################
'

FUNCTION sdxAlt ALIAS "sdxAlt" (timex AS ASCIIZ,BYVAL Longitude AS DOUBLE, BYVAL Latitude AS DOUBLE) EXPORT AS DOUBLE

IF xgetv($TagError,Timex) THEN
   FUNCTION = %BadParameter
   EXIT FUNCTION
END IF
IF Longitude = %BadParameter THEN
   FUNCTION = %BadParameter
   EXIT FUNCTION
END IF
IF Latitude = %BadParameter THEN
   FUNCTION = %BadParameter
   EXIT FUNCTION
END IF

DIM xid AS STRING
DIM xvalue AS EXT
xvalue = xUnWrapTime(xid,Timex)
IF xvalue = %BadParameter THEN
   FUNCTION = %BadParameter
   EXIT FUNCTION
END IF

IF BadLatitude(Latitude) THEN
   FUNCTION = %BadParameter
   EXIT FUNCTION
END IF
IF BadLongitude(Longitude) THEN
   FUNCTION = %BadParameter
   EXIT FUNCTION
END IF

SELECT CASE xid
       CASE $TagMJD
            IF BadMJD(xvalue) THEN
               FUNCTION = %BadParameter
               EXIT FUNCTION
            END IF
            FUNCTION = INTAltMJDUnit(xvalue, Longitude , Latitude, %Unitdegree)
            EXIT FUNCTION
       CASE $TagSDY
            IF BadSDY(xvalue) THEN
               FUNCTION = %BadParameter
               EXIT FUNCTION
            END IF
            FUNCTION = INTAltSDYUnit(xvalue, Longitude , Latitude, %Unitdegree)
            EXIT FUNCTION
       CASE ELSE
            FUNCTION = %BadParameter
            EXIT FUNCTION
END SELECT
END FUNCTION


'
' ########################
' sdFindAnnFeatx ()
' ########################
'

FUNCTION sdxFindAnnFeatx ALIAS "sdxFindAnnFeatx" (timex AS ASCIIZ, BYVAL Feature AS LONG) EXPORT AS STRING

' Check Bad Parm

DIM Scale AS LONG
DIM Year AS LONG

IF xgetv($TagError,Timex) THEN
   FUNCTION = xerror(%BadParameter)
   EXIT FUNCTION
END IF

DIM xid AS STRING
DIM xvalue AS EXT
xvalue = xUnWrapTime(xid,Timex)

IF xvalue = %BadParameter THEN
   FUNCTION = xerror(%BadParameter)
   EXIT FUNCTION
END IF

SELECT CASE xid
       CASE $TagMJD
            IF BadMJD(xvalue) THEN
               FUNCTION = xerror(%BadParameter)
               EXIT FUNCTION
            END IF
            Scale = %ScaleMJD
            PartsOfJD MJD2JD(xvalue),Year,0,0,0,0,0
        CASE $TagSDY
            IF BadSDY(xvalue) THEN
               FUNCTION = xerror(%BadParameter)
               EXIT FUNCTION
            END IF
            Scale = %ScaleSDY
            Year = 0
        CASE ELSE
           FUNCTION = xerror(%BadParameter)
           EXIT FUNCTION
END SELECT

IF BadFeature(Feature) THEN
   FUNCTION = xerror(%BadParameter)
   EXIT FUNCTION
END IF


DIM X AS DOUBLE
X = INTFindAnnFeat(Scale, Year, Feature)
DIM value AS STRING
value = STR$(X)
DIM TimeZone AS STRING
TimeZone = xgets($TagTimeZone, Timex)
DIM TimeZonex AS STRING
IF xgetv($TagError,TimeZone) THEN
    TimeZonex = ""
ELSE
    TimeZonex = xwrap($TagTimeZone,TimeZone)
END IF


SELECT CASE xid
        CASE $TagMJD
            FUNCTION = value+xwrap($TagMJD,value)+TimeZonex
        CASE $TagSDY
            FUNCTION = value+xwrap($TagSDY,value)+TimeZonex
        CASE ELSE
           FUNCTION = xerror(%BadParameter)
           EXIT FUNCTION
END SELECT

END FUNCTION

'
' ########################
' sdxTime2Text ()
' ########################
'

FUNCTION sdxTime2Text ALIAS "sdxTime2Text" (timex AS ASCIIZ) EXPORT AS STRING

IF xgetv($TagError,Timex) THEN
   FUNCTION = xerror(%BadParameter)
   EXIT FUNCTION
END IF

DIM xid AS STRING
DIM xvalue AS EXT
xvalue = xUnWrapTime(xid,Timex)
IF xvalue = %BadParameter THEN
   FUNCTION = xerror(%BadParameter)
   EXIT FUNCTION
END IF

DIM TimeZone AS DOUBLE
TimeZone = xGetv($TagTimeZone,Timex)

DIM year AS LONG
DIM month AS LONG
DIM day AS LONG
DIM hour AS LONG
DIM minute AS LONG
DIM sec AS EXT

DIM Ans AS STRING

SELECT CASE xid
   CASE $TagCTD
        IF BadCTD(xvalue) THEN
           FUNCTION = xerror(%BadParameter)
           EXIT FUNCTION
        END IF
        PartsOfJD MJD2JD(ymd2mjd(2001,1,1,0,0,0,0)+xvalue),0,0,0,Hour,Minute,Sec
        Sec = ROUND(Sec,5)
        Ans = FORMAT$(hour,"#0")+":"+FORMAT$(minute,"#0")+":"+FORMAT$(sec,"#0.####")
        FUNCTION = Ans
        EXIT FUNCTION
   CASE $TagUTD
        IF BadUTD(xvalue) THEN
           FUNCTION = xerror(%BadParameter)
           EXIT FUNCTION
        END IF
        PartsOfJD MJD2JD(ymd2mjd(2001,1,1,0,0,0,-TimeZone)+xvalue),0,0,0,Hour,Minute,Sec
        Sec = ROUND(Sec,5)
        Ans = FORMAT$(hour,"#0")+":"+FORMAT$(minute,"#0")+":"+FORMAT$(sec,"#0.####")
   CASE $TagSDY
        IF BadSDY(xvalue) THEN
           FUNCTION = xerror(%BadParameter)
           EXIT FUNCTION
        END IF
        PartsOfJD MJD2JD(SDY2MJD(2001,xvalue))-TimeZone/24,0,Month,Day,Hour,Minute,Sec
        Sec = ROUND(Sec,5)
        Ans = FORMAT$(month,"#0")+"/"+FORMAT$(day,"#0")+" "+FORMAT$(hour,"#0")+":"+FORMAT$(minute,"#0")+":"+FORMAT$(sec,"#0.####")
   CASE $TagMJD
        IF BadMJD(xvalue) THEN
           FUNCTION = xerror(%BadParameter)
           EXIT FUNCTION
        END IF
        PartsOfJD MJD2JD(xvalue)-TimeZone/24,Year,Month,Day,Hour,Minute,Sec
        Sec = ROUND(Sec,5)
        Ans = FORMAT$(year,"#####")+"/"+FORMAT$(month,"#0")+"/"+FORMAT$(day,"#0")+" "+FORMAT$(hour,"#0")+":"+FORMAT$(minute,"#0")+":"+FORMAT$(sec,"#0.####")
   CASE ELSE
        FUNCTION = xerror(%BadParameter)
        EXIT FUNCTION
END SELECT
IF TimeZone<> 0 THEN
    DIM TimeZones AS STRING
    TimeZones=FORMAT$(TimeZone,"")
    IF TimeZone > 0 THEN
        TimeZones = "+"+TimeZones
    END IF

   Ans = Ans + " ["+TimeZones+" gives UT]"
ELSE
    Ans = Ans + " [UT]"
END IF
FUNCTION = Ans

END FUNCTION

'
' ########################
' sdDUnpackx ()
' ########################
'
FUNCTION sdDUnpackx ALIAS "sdDUnpackx" (BYVAL Degree AS DOUBLE) EXPORT AS STRING

IF Degree = %BadParameter THEN
   FUNCTION = xerror(%BadParameter)
   EXIT FUNCTION
END IF

DIM Sign AS LONG
DIM D AS LONG
DIM M AS LONG
DIM S AS DOUBLE
PartsOfD Degree,Sign,D,M,S
S = ROUND(S,10)
FUNCTION = xwrap($TagDMSSign,STR$(Sign))+xwrap($TagDegree,STR$(D))+xwrap($TagArcMinute,STR$(M))+xwrap($TagArcSec,STR$(S))
END FUNCTION

'
' ########################
' sdD2Text ()
' ########################
'
FUNCTION sdD2Text ALIAS "sdD2Text" (BYVAL Degree AS DOUBLE) EXPORT AS STRING


IF Degree = %BadParameter THEN
   FUNCTION = xerror(%BadParameter)
   EXIT FUNCTION
END IF

DIM Sign AS LONG
DIM D AS LONG
DIM M AS LONG
DIM S AS DOUBLE
PartsOfD Degree,Sign,D,M,S
S = ROUND(S,4)
DIM SignSymbol$
IF sign = 1 THEN
   SignSymbol$ = " "
ELSE
   SignSymbol$ = "-"
END IF
FUNCTION = SignSymbol$+FORMAT$(D,"0")+":"+FORMAT$(M,"#0")+":"+FORMAT$(S,"#0.####")
END FUNCTION
'
' ########################
' sdxTimeUnpackx ()
' ########################
'

FUNCTION sdxTimeUnpackx ALIAS "sdxTimeUnpackx" (timex AS ASCIIZ) EXPORT AS STRING

IF xgetv($TagError,Timex) THEN
   FUNCTION = xerror(%BadParameter)
   EXIT FUNCTION
END IF

DIM xid AS STRING
DIM xvalue AS EXT
xvalue = xUnWrapTime(xid,Timex)
IF xvalue = %BadParameter THEN
   FUNCTION = xerror(%BadParameter)
   EXIT FUNCTION
END IF

DIM TimeZone AS DOUBLE
TimeZone = xGetv($TagTimeZone,Timex)

DIM year AS LONG
DIM month AS LONG
DIM day AS LONG
DIM hour AS LONG
DIM minute AS LONG
DIM sec AS EXT

DIM Ans AS STRING

SELECT CASE xid
   CASE $TagCTD
        IF BadCTD(xvalue) THEN
           FUNCTION = xerror(%BadParameter)
           EXIT FUNCTION
        END IF
        PartsOfJD MJD2JD(ymd2mjd(2001,1,1,0,0,0,0)+xvalue),0,0,0,Hour,Minute,Sec
        Sec = ROUND(Sec,5)
        Ans = xwrap($TagHour,STR$(hour))+xwrap($TagTMinute,STR$(minute))+xwrap($TagTSec,STR$(sec))
        FUNCTION = Ans
        EXIT FUNCTION
   CASE $TagUTD
        IF BadUTD(xvalue) THEN
           FUNCTION = xerror(%BadParameter)
           EXIT FUNCTION
        END IF
        PartsOfJD MJD2JD(ymd2mjd(2001,1,1,0,0,0,-TimeZone)+xvalue),0,0,0,Hour,Minute,Sec
        Sec = ROUND(Sec,5)
        Ans = xwrap($TagHour,STR$(hour))+xwrap($TagTMinute,STR$(minute))+xwrap($TagTSec,STR$(sec))
   CASE $TagSDY
        IF BadSDY(xvalue) THEN
           FUNCTION = xerror(%BadParameter)
           EXIT FUNCTION
        END IF
        PartsOfJD MJD2JD(SDY2MJD(2001,xvalue))-TimeZone/24,0,Month,Day,Hour,Minute,Sec
        Sec = ROUND(Sec,5)
        Ans = xwrap($TagMonth,STR$(Month))+xwrap($Tagday,STR$(day))+xwrap($TagHour,STR$(hour))+xwrap($TagTMinute,STR$(minute))+xwrap($TagTSec,STR$(sec))
   CASE $TagMJD
        IF BadMJD(xvalue) THEN
           FUNCTION = xerror(%BadParameter)
           EXIT FUNCTION
        END IF
        PartsOfJD MJD2JD(xvalue)-TimeZone/24,Year,Month,Day,Hour,Minute,Sec
        Sec = ROUND(Sec,5)
        Ans = xwrap($TagYear,STR$(year))+xwrap($TagMonth,STR$(Month))+xwrap($Tagday,STR$(day))+xwrap($TagHour,STR$(hour))+xwrap($TagTMinute,STR$(minute))+xwrap($TagTSec,STR$(sec))
   CASE ELSE
        FUNCTION = xerror(%BadParameter)
        EXIT FUNCTION
END SELECT
IF TimeZone<> 0 THEN
   Ans = Ans + xwrap($TagTimeZone,STR$(TimeZone))
END IF
FUNCTION = Ans

END FUNCTION

'
' ########################
' sdDMS2D ()
' ########################
'
FUNCTION sdDMS2D ALIAS "sdDMS2D" (BYVAL Sign AS LONG, BYVAL Degree AS DOUBLE,BYVAL Minute AS DOUBLE, BYVAL Sec AS DOUBLE) EXPORT AS DOUBLE

IF Sign = %BadParameter THEN
    FUNCTION = %BadParameter
    EXIT FUNCTION
END IF
IF Degree = %BadParameter THEN
    FUNCTION = %BadParameter
    EXIT FUNCTION
END IF
IF Minute = %BadParameter THEN
    FUNCTION = %BadParameter
    EXIT FUNCTION
END IF
IF Sec = %BadParameter THEN
    FUNCTION = %BadParameter
    EXIT FUNCTION
END IF

IF Sign <> 1 AND Sign <> -1 THEN
    FUNCTION = %BadParameter
    EXIT FUNCTION
END IF
IF Degree < 0 THEN
    FUNCTION = %BadParameter
    EXIT FUNCTION
END IF
IF Minute < 0 THEN
    FUNCTION = %BadParameter
    EXIT FUNCTION
END IF
IF Sec < 0 THEN
    FUNCTION = %BadParameter
    EXIT FUNCTION
END IF

FUNCTION = Sign*DMS2Unit((Degree),(Minute),(Sec),%UnitDegree)
END FUNCTION

'
' ########################
' sdD2Unit ()
' ########################
'
FUNCTION sdD2Unit ALIAS "sdD2Unit" (BYVAL Degree AS DOUBLE, BYVAL ToUnit AS LONG) EXPORT AS DOUBLE
' debug code
'Dim i&
'For i& = 1 To Len(ToUnit)
'    MsgBox Str$(Asc(Tounit,i&))
'Next
FUNCTION = INTUnit2Unit(Degree, %UnitDegree, ToUnit)
END FUNCTION
'
' ########################
' sdUnit2D ()
' ########################
'
FUNCTION sdUnit2D ALIAS "sdUnit2D" (BYVAL Value AS DOUBLE, BYVAL FromUnit AS LONG) EXPORT AS DOUBLE
FUNCTION = INTUnit2Unit(Value, FromUnit, %UnitDegree)
END FUNCTION




'
' ########################
' sdEqu2Azi ()
' ########################
'
FUNCTION sdEqu2Azi ALIAS "sdEqu2Azi" (BYVAL H AS DOUBLE,BYVAL Latitude AS DOUBLE, BYVAL  Decl AS DOUBLE) EXPORT AS DOUBLE

IF H = %BadParameter THEN
   FUNCTION = %BadParameter
   EXIT FUNCTION
END IF
IF Latitude = %BadParameter THEN
   FUNCTION = %BadParameter
   EXIT FUNCTION
END IF
IF Decl = %BadParameter THEN
   FUNCTION = %BadParameter
   EXIT FUNCTION
END IF

IF BadLatitude(Latitude) THEN
   FUNCTION = %BadParameter
   EXIT FUNCTION
END IF
IF BadDecl(Decl) THEN
   FUNCTION = %BadParameter
   EXIT FUNCTION
END IF

FUNCTION = INTHPD2Azi(H,Latitude,Decl)

END FUNCTION

'
' ########################
' sdEqu2Alt ()
' ########################
'
FUNCTION sdEqu2Alt ALIAS "sdEqu2Alt" (BYVAL H AS DOUBLE,BYVAL Latitude AS DOUBLE, BYVAL  Decl AS DOUBLE) EXPORT AS DOUBLE

IF H = %BadParameter THEN
   FUNCTION = %BadParameter
   EXIT FUNCTION
END IF
IF Latitude = %BadParameter THEN
   FUNCTION = %BadParameter
   EXIT FUNCTION
END IF
IF Decl = %BadParameter THEN
   FUNCTION = %BadParameter
   EXIT FUNCTION
END IF

IF BadLatitude(Latitude) THEN
   FUNCTION = %BadParameter
   EXIT FUNCTION
END IF
IF BadDecl(Decl) THEN
   FUNCTION = %BadParameter
   EXIT FUNCTION
END IF

FUNCTION = INTHPD2Alt(H,Latitude,Decl)

END FUNCTION

'
' ########################
' sdSin ()
' ########################
'
FUNCTION sdSin ALIAS "sdSin" (BYVAL value AS DOUBLE ) EXPORT AS DOUBLE
IF value = %BadParameter THEN
   FUNCTION = %BadParameter
   EXIT FUNCTION
END IF
FUNCTION = dSin((value))
END FUNCTION

'
' ########################
' sdCos ()
' ########################
'
FUNCTION sdCos ALIAS "sdCos" (BYVAL value AS DOUBLE ) EXPORT AS DOUBLE
IF value = %BadParameter THEN
   FUNCTION = %BadParameter
   EXIT FUNCTION
END IF
FUNCTION = dCos((value))
END FUNCTION

'
' ########################
' sdTan ()
' ########################
'
FUNCTION sdTan ALIAS "sdTan" (BYVAL value AS DOUBLE ) EXPORT AS DOUBLE
IF value = %BadParameter THEN
   FUNCTION = %BadParameter
   EXIT FUNCTION
END IF
FUNCTION = dTan((value))
END FUNCTION

'
' ########################
' sdATan ()
' ########################
'
FUNCTION sdATan ALIAS "sdATan" (BYVAL value AS DOUBLE) EXPORT AS DOUBLE
IF value = %BadParameter THEN
   FUNCTION = %BadParameter
   EXIT FUNCTION
END IF
FUNCTION = dAtan((value))
END FUNCTION

'
' ########################
' sdASin ()
' ########################
'
FUNCTION sdASin ALIAS "sdASin" (BYVAL value AS DOUBLE ) EXPORT AS DOUBLE
IF value = %BadParameter THEN
   FUNCTION = %BadParameter
   EXIT FUNCTION
END IF
IF value < -1 OR 1 < value THEN
   FUNCTION = %BadParameter
   EXIT FUNCTION
END IF
FUNCTION = dAsin((value))
END FUNCTION

'
' ########################
' sdACos ()
' ########################
'
FUNCTION sdACos ALIAS "sdACos" (BYVAL value AS DOUBLE ) EXPORT AS DOUBLE
IF value = %BadParameter THEN
   FUNCTION = %BadParameter
   EXIT FUNCTION
END IF
IF value < -1 OR 1 < value THEN
   FUNCTION = %BadParameter
   EXIT FUNCTION
END IF
FUNCTION = dAcos((value))
END FUNCTION


'
' ########################
' sdPI ()
' ########################
'
FUNCTION sdPI ALIAS "sdPI" () EXPORT AS DOUBLE
FUNCTION  = INTPI() ' 4 * Atn(1)
END FUNCTION


'
' ########################
' sdXval ()
' ########################
'
FUNCTION sdxGetv ALIAS "sdxGetv" (selector AS ASCIIZ, xstring AS ASCIIZ) EXPORT AS DOUBLE

FUNCTION = xGetv(selector, xstring)

END FUNCTION










    